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,36 +2,32 @@
 # See LICENSE.txt for details.
 
 import logging
-import warnings
 from pathlib import Path
-from typing import Optional, Union
+from typing import Optional, Sequence, Union
 
-import matplotlib.pyplot as plt
+import click
 import pandas as pd
 import yaml
 
 from spectre.Pipelines.Bbh.InitialData import generate_id
-from spectre.Pipelines.EccentricityControl.EccentricityControl import (
-    coordinate_separation_eccentricity_control,
-)
-
-# Suppress specific RuntimeWarnings
-warnings.filterwarnings(
-    "ignore", message="Number of calls to function has reached maxfev"
+from spectre.Pipelines.EccentricityControl.EccentricityControlParams import (
+    eccentricity_control_params,
+    eccentricity_control_params_options,
 )
+from spectre.support.Schedule import scheduler_options
 
 logger = logging.getLogger(__name__)
 
 
 def eccentricity_control(
-    h5_file,
-    id_input_file_path,
+    h5_files: Union[Union[str, Path], Sequence[Union[str, Path]]],
+    id_input_file_path: Union[str, Path],
     pipeline_dir: Union[str, 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
+    # Eccentricity control parameters
+    tmin: Optional[float] = 500,
+    tmax: Optional[float] = None,
+    plot_output_dir: Optional[Union[str, Path]] = None,
+    # Scheduler options
     **scheduler_kwargs,
 ):
     """Eccentricity reduction post inspiral.
@@ -47,9 +43,9 @@ def eccentricity_control(
       at 500 and using all available data by default, with the option to adjust
       'tmin' and 'tmax' dynamically.
 
-    - Calls the 'coordinate_separation_eccentricity_control' function to
-      calculate the current eccentricity and extract more accurate orbital
-      parameters.
+    - Get the new orbital parameters by calling the function
+      'eccentricity_control_params' in
+      'spectre.Pipelines.EccentricityControl.EccentricityControl'.
 
     - Displays the fit results in a tabular format using a pandas DataFrame.
 
@@ -60,23 +56,17 @@ def eccentricity_control(
       'generate_id' function.
 
     Arguments:
-      h5_file: file that contains the trajectory data
+      h5_files: files that contain the trajectory data
       id_input_file_path: path to the input file of the initial data run
       pipeline_dir : directory where the pipeline outputs are stored.
-      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
+
+    See the 'eccentricity_control_params' function for details on the other
+    arguments, as well as the 'schedule' function for the scheduling options.
     """
     # Read and process the initial data input file
     with open(id_input_file_path, "r") as open_input_file:
         id_metadata, 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"]
     id_params = id_metadata["Next"]["With"]
     control_params = id_params["control_params"]
     mass_A = control_params["mass_A"]
@@ -86,64 +76,44 @@ def eccentricity_control(
     x_B, x_A = binary_data["XCoords"]
     separation = x_A - x_B
 
-    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,
+    # Find the current eccentricity and determine new parameters to put into
+    # generate-id
+    (
+        eccentricity,
+        ecc_std_dev,
+        new_orbital_params,
+    ) = eccentricity_control_params(
+        h5_files,
+        id_input_file_path,
         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
-
-    fit_result = ecout["fit result"]
+        plot_output_dir=plot_output_dir,
+    )
 
-    # 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)
 
-    if fit_result["eccentricity"] < 0.001:
+    # Stop eccentricity control if eccentricity is below threshold
+    if eccentricity < 0.001:
         print("Success")
         # Should continue the simulation either by restarting from a
         # checkpoint, or from the volume data - will do later
@@ -157,10 +127,8 @@ def eccentricity_control(
         dimensionless_spin_b=spin_B,
         # Orbital parameters
         separation=separation,
-        orbital_angular_velocity=fit_result["updated xcts values"]["omega"],
-        radial_expansion_velocity=fit_result["updated xcts values"][
-            "expansion"
-        ],
+        orbital_angular_velocity=new_orbital_params["Omega0"],
+        radial_expansion_velocity=new_orbital_params["adot0"],
         # Scheduling options
         control=id_params["control"],
         refinement_level=id_params["control_refinement_level"],
@@ -170,3 +138,20 @@ def eccentricity_control(
         pipeline_dir=pipeline_dir,
         **scheduler_kwargs,
     )
+
+
+@click.command(name="eccentricity-control", help=eccentricity_control.__doc__)
+@eccentricity_control_params_options
+@click.option(
+    "--pipeline-dir",
+    "-d",
+    type=click.Path(
+        writable=True,
+        path_type=Path,
+    ),
+    help="Directory where steps in the pipeline are created.",
+)
+@scheduler_options
+def eccentricity_control_command(**kwargs):
+    _rich_traceback_guard = True  # Hide traceback until here
+    eccentricity_control(**kwargs)

support/Pipelines/Bbh/Inspiral.yaml

@@ -19,10 +19,8 @@ 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
+    plot_output_dir: ./
     id_input_file_path: {{id_input_file_path }}
     pipeline_dir: {{ pipeline_dir }}
     scheduler: {{ scheduler | default("None") }}

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
 )