Merge pull request #6431 from AlexCarpenter46/ringdown_cleanup

Clean up ringdown script

src/Evolution/Ringdown/Python/Bindings.cpp

@@ -30,7 +30,7 @@ void bind_strahlkorper_coefs_in_ringdown_distorted_frame(py::module& m) {
         py::arg("requested_number_of_times_from_end"), py::arg("match_time"),
         py::arg("settling_timescale"), py::arg("exp_func_and_2_derivs"),
         py::arg("exp_outer_bdry_func_and_2_derivs"),
-        py::arg("rot_func_and_2_derivs"));
+        py::arg("rot_func_and_2_derivs"), py::arg("trans_func_and_2_derivs"));
 }
 }  // namespace evolution::Ringdown::py_bindings
 

src/Evolution/Ringdown/Python/CMakeLists.txt

@@ -12,7 +12,6 @@ spectre_python_add_module(
     PYTHON_FILES
     __init__.py
     ComputeAhCCoefsInRingdownDistortedFrame.py
-    FunctionsOfTimeFromVolume.py
 )
 
 spectre_python_headers(

src/Evolution/Ringdown/Python/ComputeAhCCoefsInRingdownDistortedFrame.py

@@ -70,9 +70,7 @@ def fit_to_a_cubic(times, coefs, match_time, zero_coefs_eps):
 def compute_ahc_coefs_in_ringdown_distorted_frame(
     ahc_reductions_path,
     ahc_subfile,
-    exp_func_and_2_derivs,
-    exp_outer_bdry_func_and_2_derivs,
-    rot_func_and_2_derivs,
+    evaluated_fot_dict,
     number_of_ahc_finds_for_fit,
     match_time,
     settling_timescale,
@@ -119,9 +117,10 @@ def compute_ahc_coefs_in_ringdown_distorted_frame(
             number_of_ahc_finds_for_fit,
             match_time,
             settling_timescale,
-            exp_func_and_2_derivs,
-            exp_outer_bdry_func_and_2_derivs,
-            rot_func_and_2_derivs,
+            evaluated_fot_dict["Expansion"],
+            evaluated_fot_dict["ExpansionOuterBoundary"],
+            evaluated_fot_dict["Rotation"],
+            None,
         )
     )
 
@@ -151,24 +150,24 @@ def compute_ahc_coefs_in_ringdown_distorted_frame(
             -number_of_ahc_finds_for_fit:
         ]
 
-    logger.info("AhC times available: " + str(ahc_times.shape[0]))
-    logger.info(
+    logger.debug("AhC times available: " + str(ahc_times.shape[0]))
+    logger.debug(
         "AhC available time range: "
         + str(np.min(ahc_times))
         + " - "
         + str(np.max(ahc_times))
     )
-    logger.info("AhC times used: " + str(ahc_times_for_fit.shape[0]))
-    logger.info(
+    logger.debug("AhC times used: " + str(ahc_times_for_fit.shape[0]))
+    logger.debug(
         "AhC used time range: "
         + str(np.min(ahc_times_for_fit))
         + " - "
         + str(np.max(ahc_times_for_fit))
     )
-    logger.info(
+    logger.debug(
         "Coef times available: " + str(coefs_at_different_times.shape[0])
     )
-    logger.info(
+    logger.debug(
         "Coef times used: " + str(coefs_at_different_times_for_fit.shape[0])
     )
 

src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp

@@ -26,9 +26,13 @@ std::vector<DataVector> strahlkorper_coefs_in_ringdown_distorted_frame(
     const std::string& surface_subfile_name,
     const size_t requested_number_of_times_from_end, const double match_time,
     const double settling_timescale,
-    const std::array<double, 3>& exp_func_and_2_derivs,
-    const std::array<double, 3>& exp_outer_bdry_func_and_2_derivs,
-    const std::vector<std::array<double, 4>>& rot_func_and_2_derivs) {
+    const std::optional<std::array<double, 3>>& exp_func_and_2_derivs,
+    const std::optional<std::array<double, 3>>&
+        exp_outer_bdry_func_and_2_derivs,
+    const std::optional<std::vector<std::array<double, 4>>>&
+        rot_func_and_2_derivs,
+    const std::optional<std::array<std::array<double, 3>, 3>>&
+        trans_func_and_2_derivs) {
   // Read the AhC coefficients from the H5 file
   const std::vector<ylm::Strahlkorper<Frame::Inertial>>& ahc_inertial_h5 =
       ylm::read_surface_ylm<Frame::Inertial>(
@@ -49,16 +53,35 @@ std::vector<DataVector> strahlkorper_coefs_in_ringdown_distorted_frame(
   // Create a time-dependent domain; only the the time-dependent map options
   // matter; the domain is just a spherical shell with inner and outer
   // radii chosen so any conceivable common horizon will fit between them.
-  const domain::creators::time_dependent_options::ExpansionMapOptions<true>
-      expansion_map_options{exp_func_and_2_derivs, settling_timescale,
-                            exp_outer_bdry_func_and_2_derivs,
-                            settling_timescale};
-  const domain::creators::time_dependent_options::RotationMapOptions<true>
-      rotation_map_options{rot_func_and_2_derivs, settling_timescale};
+  const auto expansion_map_options =
+      exp_func_and_2_derivs.has_value()
+          ? domain::creators::time_dependent_options::
+                ExpansionMapOptions<true>{exp_func_and_2_derivs.value(),
+                                    settling_timescale,
+                                    exp_outer_bdry_func_and_2_derivs.value(),
+                                    settling_timescale}
+          : std::optional<domain::creators::time_dependent_options::
+                              ExpansionMapOptions<true>>{};
+  const auto rotation_map_options =
+      rot_func_and_2_derivs.has_value()
+          ? domain::creators::time_dependent_options::
+                RotationMapOptions<true>{rot_func_and_2_derivs.value(),
+                                   settling_timescale}
+          : std::optional<domain::creators::time_dependent_options::
+                              RotationMapOptions<true>>{};
+  const auto translation_map_options =
+      trans_func_and_2_derivs.has_value()
+          ? domain::creators::sphere::TimeDependentMapOptions::
+                TranslationMapOptions{trans_func_and_2_derivs.value()}
+          : std::optional<domain::creators::sphere::TimeDependentMapOptions::
+                              TranslationMapOptions>{};
   const domain::creators::sphere::TimeDependentMapOptions
-      time_dependent_map_options{match_time,           std::nullopt,
-                                 rotation_map_options, expansion_map_options,
-                                 std::nullopt,         true};
+      time_dependent_map_options{match_time,
+                                 std::nullopt,
+                                 rotation_map_options,
+                                 expansion_map_options,
+                                 translation_map_options,
+                                 true};
   const domain::creators::Sphere domain_creator{
       0.01,
       200.0,

src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp

@@ -37,7 +37,12 @@ std::vector<DataVector> strahlkorper_coefs_in_ringdown_distorted_frame(
     const std::string& surface_subfile_name,
     size_t requested_number_of_times_from_end, double match_time,
     double settling_timescale,
-    const std::array<double, 3>& exp_func_and_2_derivs,
-    const std::array<double, 3>& exp_outer_bdry_func_and_2_derivs,
-    const std::vector<std::array<double, 4>>& rot_func_and_2_derivs);
+    const std::optional<std::array<double, 3>>& exp_func_and_2_derivs =
+        std::nullopt,
+    const std::optional<std::array<double, 3>>&
+        exp_outer_bdry_func_and_2_derivs = std::nullopt,
+    const std::optional<std::vector<std::array<double, 4>>>&
+        rot_func_and_2_derivs = std::nullopt,
+    const std::optional<std::array<std::array<double, 3>, 3>>&
+        trans_func_and_2_derivs = std::nullopt);
 }  // namespace evolution::Ringdown

src/IO/H5/Python/CMakeLists.txt

@@ -22,6 +22,7 @@ spectre_python_add_module(
   ExtendConnectivityData.py
   ExtractDatFromH5.py
   ExtractInputSourceYamlFromH5.py
+  FunctionsOfTimeFromVolume.py
   IterElements.py
   MODULE_PATH "IO"
   )

src/IO/H5/Python/FunctionsOfTimeFromVolume.py

@@ -0,0 +1,64 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import logging
+from pathlib import Path
+from typing import Optional, Union
+
+import click
+import numpy as np
+import yaml
+from rich.pretty import pretty_repr
+
+import spectre.IO.H5 as spectre_h5
+from spectre.Domain import deserialize_functions_of_time
+
+logger = logging.getLogger(__name__)
+
+
+def functions_of_time_from_volume(
+    fot_vol_h5_path, fot_vol_subfile, match_time, fot_to_observe=None
+):
+    """This function returns a dictionary of the FunctionsOfTime from a volume
+    h5 file evaluated at a desired time along with the time closest to the
+    match_time passed in.
+
+    Arguments:
+    fot_vol_h5_path: The full path to volume data containing functions of time
+    that have been evaluated at the match time.
+    fot_vol_subfile: The subfile containing volume data with functions of time
+    evaulated at match time.
+    match_time: The desired time to retrieve the function of time values.
+    """
+
+    functions_of_time_at_match_time_dict = {}
+
+    with spectre_h5.H5File(fot_vol_h5_path, "r") as h5file:
+        if fot_vol_subfile.split(".")[-1] == "vol":
+            fot_vol_subfile = fot_vol_subfile.split(".")[0]
+        volfile = h5file.get_vol("/" + fot_vol_subfile)
+        obs_ids = volfile.list_observation_ids()
+        fot_times = np.array(list(map(volfile.get_observation_value, obs_ids)))
+        which_obs_id = np.argmin(np.abs(fot_times - match_time))
+        serialized_fots = volfile.get_functions_of_time(obs_ids[which_obs_id])
+        functions_of_time = deserialize_functions_of_time(serialized_fots)
+        logger.debug("Desired match time: " + str(match_time))
+        logger.debug("Selected ObservationID: " + str(which_obs_id))
+        logger.debug("Selected match time: " + str(fot_times[which_obs_id]))
+
+        functions_of_time_at_match_time_dict["MatchTime"] = fot_times[
+            which_obs_id
+        ]
+
+        for fot_name, fot in functions_of_time.items():
+            fot_at_match_time = fot.func_and_2_derivs(fot_times[which_obs_id])
+            if len(fot_at_match_time[0]) != 1:
+                functions_of_time_at_match_time_dict[fot_name] = [
+                    [coef for coef in x] for x in fot_at_match_time
+                ]
+            else:
+                functions_of_time_at_match_time_dict[fot_name] = [
+                    x[0] for x in fot_at_match_time
+                ]
+
+    return functions_of_time_at_match_time_dict