Add functionality to Ylm IO

src/Domain/Creators/Sphere.cpp

@@ -159,10 +159,9 @@ Sphere::Sphere(
       tnsr::I<double, 3, Frame::Grid>{std::array{0.0, 0.0, 0.0}};
 
   // Determine number of blocks
-  num_blocks_per_shell_ =
-      which_wedges_ == ShellWedges::All
-          ? 6
-          : which_wedges_ == ShellWedges::FourOnEquator ? 4 : 1;
+  num_blocks_per_shell_ = which_wedges_ == ShellWedges::All             ? 6
+                          : which_wedges_ == ShellWedges::FourOnEquator ? 4
+                                                                        : 1;
   num_blocks_ = num_blocks_per_shell_ * num_shells_ + (fill_interior_ ? 1 : 0);
 
   // Create block names and groups
@@ -310,7 +309,8 @@ Domain<3> Sphere::create_domain() const {
           inner_radius_, outer_radius_,
           fill_interior_ ? std::get<InnerCube>(interior_).sphericity : 1.0, 1.0,
           use_equiangular_map_, false, radial_partitioning_,
-          radial_distribution_, which_wedges_), compression);
+          radial_distribution_, which_wedges_),
+      compression);
 
   std::unordered_map<std::string, ExcisionSphere<3>> excision_spheres{};
 

src/Domain/Creators/SphereTimeDependentMaps.hpp

@@ -169,8 +169,7 @@ struct TimeDependentMapOptions {
     using options = tmpl::list<InitialValues, DecayTimescaleRotation>;
 
     std::array<std::array<double, 4>, 3> initial_values{};
-    double decay_timescale{
-        std::numeric_limits<double>::signaling_NaN()};
+    double decay_timescale{std::numeric_limits<double>::signaling_NaN()};
   };
 
   struct ExpansionMapOptions {

src/NumericalAlgorithms/SphericalHarmonics/IO/CMakeLists.txt

@@ -8,13 +8,15 @@ add_spectre_library(${LIBRARY})
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
+  FillYlmLegendAndData.cpp
   ReadSurfaceYlm.cpp
   )
 
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
+  FillYlmLegendAndData.hpp
   ReadSurfaceYlm.hpp
   )
 
@@ -24,7 +26,7 @@ target_link_libraries(
   DataStructures
   ErrorHandling
   H5
-  Utilities
   PUBLIC
   SphericalHarmonics
+  Utilities
   )

src/ParallelAlgorithms/Interpolation/Callbacks/ObserveSurfaceData.cpp → src/NumericalAlgorithms/SphericalHarmonics/IO/FillYlmLegendAndData.cpp

@@ -1,7 +1,7 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
-#include "ParallelAlgorithms/Interpolation/Callbacks/ObserveSurfaceData.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/IO/FillYlmLegendAndData.hpp"
 
 #include <array>
 #include <cstddef>
@@ -19,7 +19,7 @@
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeString.hpp"
 
-namespace intrp::callbacks::detail {
+namespace ylm {
 template <typename Frame>
 void fill_ylm_legend_and_data(
     const gsl::not_null<std::vector<std::string>*> legend,
@@ -75,12 +75,12 @@ void fill_ylm_legend_and_data(
     }
   }
 }
-}  // namespace intrp::callbacks::detail
+}  // namespace ylm
 
 #define FRAMETYPE(instantiation_data) BOOST_PP_TUPLE_ELEM(0, instantiation_data)
 
 #define INSTANTIATE(_, instantiation_data)                                  \
-  template void intrp::callbacks::detail::fill_ylm_legend_and_data(         \
+  template void ylm::fill_ylm_legend_and_data(                              \
       gsl::not_null<std::vector<std::string>*> legend,                      \
       gsl::not_null<std::vector<double>*> data,                             \
       const ylm::Strahlkorper<FRAMETYPE(instantiation_data)>& strahlkorper, \

src/NumericalAlgorithms/SphericalHarmonics/IO/FillYlmLegendAndData.hpp

@@ -0,0 +1,29 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+#include <string>
+#include <vector>
+
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace ylm {
+/*!
+ * \brief Fills the legend and row of spherical harmonic data to write to disk
+ *
+ * The number of coefficients to write is based on `max_l`, the maximum value
+ * that the input `strahlkorper` could possibly have. When
+ * `strahlkorper.l_max() < max_l`, coefficients with \f$l\f$ higher than
+ * `strahlkorper.l_max()` will simply be zero. Assuming the same `max_l` is
+ * always used for a given surface, we will always write the same number of
+ * columns for each row, as `max_l` sets the number of columns to write
+ */
+template <typename Frame>
+void fill_ylm_legend_and_data(gsl::not_null<std::vector<std::string>*> legend,
+                              gsl::not_null<std::vector<double>*> data,
+                              const ylm::Strahlkorper<Frame>& strahlkorper,
+                              double time, size_t max_l);
+}  // namespace ylm

src/NumericalAlgorithms/SphericalHarmonics/IO/ReadSurfaceYlm.cpp

@@ -6,6 +6,7 @@
 #include <array>
 #include <cmath>
 #include <cstddef>
+#include <cstdint>
 #include <numeric>
 #include <string>
 #include <vector>
@@ -18,6 +19,7 @@
 #include "IO/H5/File.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/SpherepackIterator.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+#include "Utilities/EqualWithinRoundoff.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
@@ -152,6 +154,51 @@ Strahlkorper<Frame> read_surface_ylm_row(const Matrix& ylm_data,
 }
 }  // namespace
 
+template <typename Frame>
+ylm::Strahlkorper<Frame> read_surface_ylm_single_time(
+    const std::string& file_name, const std::string& surface_subfile_name,
+    const double time, const double relative_epsilon) {
+  h5::H5File<h5::AccessType::ReadOnly> file{file_name};
+  const std::string ylm_subfile_name{std::string{"/"} + surface_subfile_name};
+  const auto& ylm_file = file.get<h5::Dat>(ylm_subfile_name);
+  const auto& ylm_legend = ylm_file.get_legend();
+  check_legend<Frame>(ylm_legend);
+
+  std::vector<size_t> columns(ylm_legend.size());
+  std::iota(std::begin(columns), std::end(columns), 0);
+
+  const auto& dimensions = ylm_file.get_dimensions();
+  const auto num_rows = static_cast<size_t>(dimensions[0]);
+  const Matrix times =
+      ylm_file.get_data_subset(std::vector<size_t>{0_st}, 0, num_rows);
+
+  std::optional<size_t> row_number{};
+  for (size_t i = 0; i < num_rows; i++) {
+    if (equal_within_roundoff(time, times(i, 0), relative_epsilon, time)) {
+      if (row_number.has_value()) {
+        ERROR("Found more than one time in the subfile "
+              << surface_subfile_name << " of the H5 file " << file_name
+              << " that is within a relative_epsilon of " << relative_epsilon
+              << " of the time requested " << time);
+      }
+
+      row_number = i;
+    }
+  }
+
+  if (not row_number.has_value()) {
+    ERROR("Could not find time " << time << " in subfile "
+                                 << surface_subfile_name << " of H5 file "
+                                 << file_name << ". Available times are:\n"
+                                 << times);
+  }
+
+  const Matrix data = ylm_file.get_data_subset(columns, row_number.value());
+
+  // Zero for row number because this matrix only has one row
+  return read_surface_ylm_row<Frame>(data, 0);
+}
+
 template <typename Frame>
 std::vector<Strahlkorper<Frame>> read_surface_ylm(
     const std::string& file_name, const std::string& surface_subfile_name,
@@ -200,11 +247,15 @@ std::vector<Strahlkorper<Frame>> read_surface_ylm(
 
 #define FRAMETYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
 
-#define INSTANTIATE(_, data)                                       \
-  template std::vector<ylm::Strahlkorper<FRAMETYPE(data)>>         \
-  ylm::read_surface_ylm<>(const std::string& file_name,            \
-                          const std::string& surface_subfile_name, \
-                          size_t requested_number_of_times_from_end);
+#define INSTANTIATE(_, data)                                                   \
+  template std::vector<ylm::Strahlkorper<FRAMETYPE(data)>>                     \
+  ylm::read_surface_ylm<>(const std::string& file_name,                        \
+                          const std::string& surface_subfile_name,             \
+                          size_t requested_number_of_times_from_end);          \
+  template ylm::Strahlkorper<FRAMETYPE(data)>                                  \
+  ylm::read_surface_ylm_single_time<>(const std::string& file_name,            \
+                                      const std::string& surface_subfile_name, \
+                                      double time, double relative_epsilon);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE,
                         (Frame::Grid, Frame::Inertial, Frame::Distorted))

src/NumericalAlgorithms/SphericalHarmonics/IO/ReadSurfaceYlm.hpp

@@ -30,4 +30,29 @@ template <typename Frame>
 std::vector<ylm::Strahlkorper<Frame>> read_surface_ylm(
     const std::string& file_name, const std::string& surface_subfile_name,
     size_t requested_number_of_times_from_end);
+
+/*!
+ * \brief Similar to `ylm::read_surface_ylm`, this reads in spherical harmonic
+ * data for a surface and constructs a `ylm::Strahlkorper`. However, this
+ * function only does it at a specific time and returns a single
+ * `ylm::Strahlkorper`.
+ *
+ * \note If two times are found within \p epsilon of the \p time, then an error
+ * will occur. Similarly, if no \p time is found within the \p epsilon, then an
+ * error will occur as well.
+ *
+ * \param file_name name of the H5 file containing the surface's spherical
+ * harmonic data
+ * \param surface_subfile_name name of the subfile (with no leading slash nor
+ * the `.dat` extension) within `file_name` that contains the surface's
+ * spherical harmonic data to read in
+ * \param time Time to read the coefficients at.
+ * \param relative_epsilon How much error is allowed when looking for a specific
+ * time. This is useful so users don't have to know the specific time to machine
+ * precision.
+ */
+template <typename Frame>
+ylm::Strahlkorper<Frame> read_surface_ylm_single_time(
+    const std::string& file_name, const std::string& surface_subfile_name,
+    double time, double relative_epsilon);
 }  // namespace ylm

src/ParallelAlgorithms/Interpolation/CMakeLists.txt

@@ -43,6 +43,7 @@ target_link_libraries(
   Serialization
   Spectral
   SphericalHarmonics
+  SphericalHarmonicsIO
   Utilities
   INTERFACE
   Actions

src/ParallelAlgorithms/Interpolation/Callbacks/CMakeLists.txt

@@ -1,12 +1,6 @@
 # Distributed under the MIT License.
 # See LICENSE.txt for details.
 
-spectre_target_sources(
-  ${LIBRARY}
-  PRIVATE
-  ObserveSurfaceData.cpp
-  )
-
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src

src/ParallelAlgorithms/Interpolation/Callbacks/ObserveSurfaceData.hpp

@@ -17,6 +17,7 @@
 #include "IO/Observer/VolumeActions.hpp"
 #include "NumericalAlgorithms/Spectral/Basis.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/IO/FillYlmLegendAndData.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
@@ -34,22 +35,6 @@
 
 namespace intrp {
 namespace callbacks {
-namespace detail {
-// Fills the legend and row of spherical harmonic data to write to disk
-//
-// The number of coefficients to write is based on `max_l`, the maximum value
-// that the input `strahlkorper` could possibly have. When
-// `strahlkorper.l_max() < max_l`, coefficients with \f$l\f$ higher than
-// `strahlkorper.l_max()` will simply be zero. Assuming the same `max_l` is
-// always used for a given surface, we will always write the same number of
-// columns for each row, as `max_l` sets the number of columns to write
-template <typename Frame>
-void fill_ylm_legend_and_data(gsl::not_null<std::vector<std::string>*> legend,
-                              gsl::not_null<std::vector<double>*> data,
-                              const ylm::Strahlkorper<Frame>& strahlkorper,
-                              double time, size_t max_l);
-}  // namespace detail
-
 /// \brief post_interpolation_callback that outputs 2D "volume" data on a
 /// surface and the surface's spherical harmonic data
 ///
@@ -173,9 +158,9 @@ struct ObserveSurfaceData
     // coefficients regardless of the current value of l_max and (b) write a
     // constant number of columns for each row of data regardless of the current
     // l_max.
-    detail::fill_ylm_legend_and_data(make_not_null(&ylm_legend),
-                                     make_not_null(&ylm_data), strahlkorper,
-                                     time, strahlkorper.l_max());
+    ylm::fill_ylm_legend_and_data(make_not_null(&ylm_legend),
+                                  make_not_null(&ylm_data), strahlkorper, time,
+                                  strahlkorper.l_max());
 
     const std::string ylm_subfile_name{std::string{"/"} + surface_name +
                                        "_Ylm"};

tests/Unit/NumericalAlgorithms/SphericalHarmonics/IO/CMakeLists.txt

@@ -4,6 +4,7 @@
 set(LIBRARY "Test_SphericalHarmonicsIO")
 
 set(LIBRARY_SOURCES
+  Test_FillYlmLegendAndData.cpp
   Test_ReadSurfaceYlm.cpp
   )