Use SpEC's eccentricity control

cmake/SetupSpec.cmake

@@ -6,7 +6,8 @@ find_package(SpEC)
 # Make SpEC scripts available in Python. These can be used until we have ported
 # them to SpECTRE.
 if (SPEC_ROOT)
-  set(PYTHONPATH "${SPEC_ROOT}/Support/Python:${PYTHONPATH}")
+  set(PYTHONPATH "${SPEC_ROOT}/Support/Python:\
+${SPEC_ROOT}/Support/DatDataManip:${PYTHONPATH}")
 endif()
 
 if (NOT SpEC_FOUND)

support/Pipelines/Bbh/EccentricityControl.py

@@ -2,118 +2,71 @@
 # See LICENSE.txt for details.
 
 import logging
-import warnings
 
-import matplotlib.pyplot as plt
+import click
 import pandas as pd
-import yaml
 
-from spectre.Pipelines.EccentricityControl.EccentricityControl import (
-    coordinate_separation_eccentricity_control,
-)
-
-# Suppress specific RuntimeWarnings
-warnings.filterwarnings(
-    "ignore", message="Number of calls to functionhas reached maxfev"
+from spectre.Pipelines.EccentricityControl.EccentricityControlParams import (
+    eccentricity_control_params,
+    eccentricity_control_params_options,
 )
 
 logger = logging.getLogger(__name__)
 
 
-def eccentricity_control(
-    h5_file,
-    id_input_file_path,
-    subfile_name_aha="ApparentHorizons/ControlSystemAhA_Centers.dat",
-    subfile_name_ahb="ApparentHorizons/ControlSystemAhB_Centers.dat",
-    tmin=500,
-    tmax=None,  # use all available data
-    output=None,  # reads from Inspiral.yaml
-):
+def eccentricity_control(h5_files, id_input_file_path, **kwargs):
     """Eccentricity reduction post inspiral.
 
     This function can be called after the inspiral has run (see the 'Next'
-    section of the Inspiral.yam file).
+    section of the Inspiral.yaml file).
 
     This function does the following:
 
-    - Reads orbital parameters from the 'id_input_file_path'.
-
-    - For now, the time boundaries are set to start at 500, and at the end of
-      the simulation to use all available data, but we will make tmin and tmax
-      more dynamic later, so the user can change those variables according to
-      interest.
+    - Get the new orbital parameters by calling the function
+      'eccentricity_control_params' in
+      'spectre.Pipelines.EccentricityControl.EccentricityControl'.
+    - Print the new orbital parameters in a tabular format. This will be updated
+      to start a new inspiral with the updated parameters.
 
     Arguments:
       h5_file: file that contains the trajectory data
       id_input_file_path: path to the input file of the initial data run
-      subfile_name_aha: subfile for black hole A; optional parameter
-      subfile_name_ahb: subfile for black hole B; optional parameter
-      tmin: starting point for eccentricity reduction script; defaults to
-        500 if not specified
-      tmax: stopping point for eccentricity reduction script; set to use all
-        available data by default
-      output: outputs to terminal plus makes pdf file, if specified
+      **kwargs: additional arguments to be forwarded to
+        'eccentricity_control_params'.
     """
-    # Read and process the initial data input file
-    with open(id_input_file_path, "r") as open_input_file:
-        _, id_input_file = yaml.safe_load_all(open_input_file)
-    binary_data = id_input_file["Background"]["Binary"]
-    orbital_angular_velocity = binary_data["AngularVelocity"]
-    radial_expansion_velocity = binary_data["Expansion"]
-
-    if output:
-        fig = plt.figure()
-    else:
-        fig = None
-
     # Find the current eccentricity and extract new parameters to put into
-    # generate-id. Call the coordinate_separation_eccentricity_control
-    # function; extract only the ["H4"] , i.e. nonlinear fit.
-    ecout = coordinate_separation_eccentricity_control(
-        h5_file=h5_file,
-        subfile_name_aha=subfile_name_aha,
-        subfile_name_ahb=subfile_name_ahb,
-        tmin=tmin,
-        tmax=tmax,
-        angular_velocity_from_xcts=orbital_angular_velocity,
-        expansion_from_xcts=radial_expansion_velocity,
-        fig=fig,
-    )["H4"]
-
-    if output:
-        fig.savefig(output)
-
-    # Print the fit result
-    # extract new parameters to put into generate_id
+    # generate-id
+    (
+        eccentricity,
+        ecc_std_dev,
+        new_orbital_params,
+    ) = eccentricity_control_params(h5_files, id_input_file_path, **kwargs)
 
-    fit_result = ecout["fit result"]
-
-    # Prepare data from fit result
+    # Create DataFrame to display data in tabular format
     data = {
         "Attribute": [
             "Eccentricity",
-            "Omega Update",
-            "Expansion Update",
-            "Previous Omega",
-            "Previous Expansion",
-            "Updated Omega",
-            "Updated Expansion",
+            "Eccentricity error",
+            "Updated Omega0",
+            "Updated adot0",
         ],
         "Value": [
-            fit_result["eccentricity"],
-            fit_result["xcts updates"]["omega update"],
-            fit_result["xcts updates"]["expansion update"],
-            fit_result["previous xcts values"]["omega"],
-            fit_result["previous xcts values"]["expansion"],
-            fit_result["updated xcts values"]["omega"],
-            fit_result["updated xcts values"]["expansion"],
+            eccentricity,
+            ecc_std_dev,
+            new_orbital_params["Omega0"],
+            new_orbital_params["adot0"],
         ],
     }
-
-    # Create DataFrame to display data in tabular format
     df = pd.DataFrame(data)
     # Print header line
     print("=" * 40)
     # Display table
     print(df.to_string(index=False))
     print("=" * 40)
+
+
+@click.command(name="eccentricity-control", help=eccentricity_control.__doc__)
+@eccentricity_control_params_options
+def eccentricity_control_command(**kwargs):
+    _rich_traceback_guard = True  # Hide traceback until here
+    eccentricity_control(**kwargs)

support/Pipelines/Bbh/Inspiral.yaml

@@ -19,11 +19,14 @@ Next:
 Next:
   Run: spectre.Pipelines.Bbh.EccentricityControl:eccentricity_control
   With:
-    h5_file: ./BbhReductions.h5
-    subfile_name_aha: ApparentHorizons/ControlSystemAhA_Centers.dat
-    subfile_name_ahb: ApparentHorizons/ControlSystemAhB_Centers.dat
-    output: ecc_control.pdf
+    h5_files: ./BbhReductions.h5
     id_input_file_path: {{id_input_file_path }}
+    subfile_name_aha_trajectories: ApparentHorizons/ControlSystemAhA_Centers.dat
+    subfile_name_ahb_trajectories: ApparentHorizons/ControlSystemAhB_Centers.dat
+    subfile_name_aha_quantities: ObservationAhA.dat
+    subfile_name_ahb_quantities: ObservationAhB.dat
+    target_eccentricity: 0.0
+    plot_output_dir: ./
 {% endif %}
 
 ---

support/Pipelines/Bbh/__init__.py

@@ -11,6 +11,7 @@
 class Bbh(click.MultiCommand):
     def list_commands(self, ctx):
         return [
+            "eccentricity-control",
             "find-horizon",
             "generate-id",
             "postprocess-id",
@@ -19,15 +20,19 @@ def list_commands(self, ctx):
         ]
 
     def get_command(self, ctx, name):
-        if name == "find-horizon":
+        if name in ["eccentricity-control", "ecc-control"]:
+            from .EccentricityControl import eccentricity_control_command
+
+            return eccentricity_control_command
+        elif name == "find-horizon":
             from .FindHorizon import find_horizon_command
 
             return find_horizon_command
-        if name == "generate-id":
+        elif name == "generate-id":
             from .InitialData import generate_id_command
 
             return generate_id_command
-        if name == "postprocess-id":
+        elif name == "postprocess-id":
             from .PostprocessId import postprocess_id_command
 
             return postprocess_id_command

support/Pipelines/EccentricityControl/CMakeLists.txt

@@ -6,6 +6,6 @@ spectre_python_add_module(
   MODULE_PATH Pipelines
   PYTHON_FILES
   __init__.py
-  EccentricityControl.py
+  EccentricityControlParams.py
   InitialOrbitalParameters.py
 )