Add files via upload

support/Pipelines/EccentricityControl/SpecOmegaDotEccRemoval.py

@@ -0,0 +1,270 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import argparse
+import os
+import sys
+import traceback
+
+import h5py
+import numpy as np
+
+
+def inertialCoordinates(file_path, dat_file_path):
+    # Reads t, inertial_coord x, intertial_coord y, intertial_coord z
+    # and returns them as a np array.
+
+    hf = h5py.File(file_path, "r")
+    controlSystem = hf[dat_file_path][:]
+
+    columns = [1, 2, 3, 4]
+    controlSystem = controlSystem[:, [col - 1 for col in columns]]
+
+    return controlSystem
+
+
+def christodoulouMass(Dir, dat_file):
+    # returns the christodoulou masses over time in the form
+    # [[t1, CM1], [t2, CM2] ...]
+    reductionData = h5py.File(
+        os.path.join(Dir, "GhBinaryBlackHoleReductionData.h5"), "r"
+    )
+    observation = reductionData[dat_file]
+    output_array = np.zeros((observation.shape[0], 2), dtype=float)
+
+    for row, col in enumerate(observation):
+        output_array[row, 0] = col[0]
+        output_array[row, 1] = col[5]
+
+    return output_array
+
+
+def dimensionfullSpins(Dir, dat_file):
+    # uses dimensionfullSpinMag to fill out an array of the form
+    # [[t1, Sx1, Sy1, Sz1], [t2, Sx2, Sy2, Sz2]...].
+    # problem: some BBH runs only have 1 entry in ObservationAhX.dat
+    # which causes errors while interpolating
+    # to 'fix' this one could duplicate entries (10 rows is usually enough)
+
+    reductionData = h5py.File(
+        os.path.join(Dir, "GhBinaryBlackHoleReductionData.h5"), "r"
+    )
+    observation = reductionData[dat_file]
+    output_array = np.zeros((len(observation), 4), dtype=float)
+
+    for row, col in enumerate(observation):
+        output_array[row, 0] = col[0]
+        output_array[row, 2] = (
+            float(col[6]) * float(col[5]) ** 2
+        )  # dimensionlessSpinMag * CM^2 = dimensionfullSpin
+        # Sx and Sy are assumed to be 0
+    return output_array
+
+
+def editHorizons(opts):
+    try:
+        reduction_path = os.path.join(
+            opts.d, "GhBinaryBlackHoleReductionData.h5"
+        )
+        reductionData = h5py.File(reduction_path, "r")
+    except Exception as e:
+        print("Error while opening 'GhBinaryBlackHoleReductionData.h5': ", e)
+        sys.exit(1)
+
+    try:
+        horizons = h5py.File(os.path.join(opts.d, "Horizons.h5"), "a")
+    except Exception as e:
+        reductionData.close()
+        print("Error while opening 'Horizons.h5': ", e)
+        sys.exit(1)
+
+    try:
+        # for GetTrajectories()
+        horizonCentersA = inertialCoordinates(
+            reduction_path, "ApparentHorizons/ControlSystemAhA_Centers.dat"
+        )
+        horizonCentersB = inertialCoordinates(
+            reduction_path, "ApparentHorizons/ControlSystemAhB_Centers.dat"
+        )
+
+        del horizons["/AhA.dir/CoordCenterInertial.dat"]
+        del horizons["/AhB.dir/CoordCenterInertial.dat"]
+        horizons.create_dataset(
+            "/AhA.dir/CoordCenterInertial.dat", data=horizonCentersA
+        )
+        horizons.create_dataset(
+            "/AhB.dir/CoordCenterInertial.dat", data=horizonCentersB
+        )
+
+        # for GetRelaxedMasses()
+        horizonMA = christodoulouMass(opts.d, "ObservationAhA.dat")
+        horizonMB = christodoulouMass(opts.d, "ObservationAhB.dat")
+
+        del horizons["AhA.dir/ChristodoulouMass.dat"]
+        del horizons["AhB.dir/ChristodoulouMass.dat"]
+        horizons.create_dataset("AhA.dir/ChristodoulouMass.dat", data=horizonMA)
+        horizons.create_dataset("AhB.dir/ChristodoulouMass.dat", data=horizonMB)
+
+        # for GetVarsFromSpinData()
+        sA = dimensionfullSpins(opts.d, "ObservationAhA.dat")
+        sB = dimensionfullSpins(opts.d, "ObservationAhB.dat")
+
+        del horizons["AhA.dir/DimensionfulInertialSpin.dat"]
+        del horizons["AhB.dir/DimensionfulInertialSpin.dat"]
+        horizons.create_dataset("AhA.dir/DimensionfulInertialSpin.dat", data=sA)
+        horizons.create_dataset("AhB.dir/DimensionfulInertialSpin.dat", data=sB)
+
+        # param.perl file that will be given to OmegaDotEccRemoval.py
+        data = {
+            "$ID_Omega0": opts.o,
+            "$ID_adot0": opts.a,
+            "$ID_d": opts.s,
+            "$ID_MA": opts.MA,
+            "$ID_MB": opts.MB,
+        }
+
+        params = open(opts.idperl, "w")
+        for key, value in data.items():
+            params.write("{} = {};\n".format(key, value))
+        params.close()
+    except Exception as e:
+        reductionData.close()
+        horizons.close()
+        print("Error while editting 'Horizons.h5': ", e)
+        traceback.print_exc()
+        sys.exit(1)
+
+    reductionData.close()
+    horizons.close()
+    return
+
+
+if __name__ == "__main__":
+    p = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawTextHelpFormatter
+    )
+    required = p.add_argument_group("required arguments")
+
+    required.add_argument(
+        "-d",
+        type=str,
+        required=True,
+        metavar="DIR",
+        help=(
+            "Directory containing 'Horizons.h5' and"
+            " 'GhBinaryBlackHoleReductionData.h5')"
+        ),
+    )
+
+    required.add_argument(
+        "-o", type=float, required=True, metavar="FLOAT", help="Omega0"
+    )
+    required.add_argument(
+        "-a", type=float, required=True, metavar="FLOAT", help="adot0"
+    )
+    required.add_argument(
+        "-s",
+        type=float,
+        required=True,
+        metavar="FLOAT",
+        help="initial_seperation",
+    )
+
+    required.add_argument(
+        "--idperl",
+        type=str,
+        metavar="DIR",
+        help="""Full path of .perl file with parameters in the proper format.
+        Must point to the file if specified.
+        Default is %cwd%ID_Params.perl')""",
+    )
+    p.add_argument(
+        "--tmin",
+        type=float,
+        metavar="FLOAT",
+        help=(
+            "Fit points with t>tmin. Default determined by FindTmin estimation"
+            " scheme."
+        ),
+    )
+    p.add_argument(
+        "--tmax",
+        type=float,
+        metavar="FLOAT",
+        help=(
+            "Fit points with t<tmax. Defaults to min(t.back(),"
+            " 5pi/Omega0+tmin)."
+        ),
+    )
+    p.add_argument(
+        "--improved_Omega0_update",
+        action="store_true",
+        help="""If specified, multiply Omega0-update by factor Omega0/omega_r.
+                   Empirically, this is recommended to improve convergence.
+                   See ticket #305 for details.""",
+    )
+    p.add_argument(
+        "--Omega500",
+        type=float,
+        metavar="FLOAT",
+        help="""If specified, compute updating formulae, such that:
+        (1) e=0 (2) Omega(t=500)=Omega500.
+        This will give changes to Omega0, adot0 and D0,
+        and all three must be applied.
+        The updating formulae are only computed for F2cos2.""",
+    )
+    p.add_argument(
+        "--no_check",
+        action="store_true",
+        help="""If specified, do not check for large or negative
+        periastron advances. This will prevent the eccentricity from
+        being set to 9.9999""",
+    )
+    p.add_argument(
+        "--no_plot",
+        action="store_true",
+        help="If specified, do not generate plots.",
+    )
+    p.add_argument(
+        "--unbounded_fits", action="store_true", help="Use old unbounded fits."
+    )
+    p.add_argument(
+        "--varpro",
+        action="store_true",
+        help="If specified, use variable projection fitting algorithm.",
+    )
+    required.add_argument(
+        "-MA", type=float, metavar="FLOAT", help="Mass of BBH A. Defualt is 0.5"
+    )
+    required.add_argument(
+        "-MB", type=float, metavar="FLOAT", help="Mass of BBH B. Default is 0.5"
+    )
+
+    opts = p.parse_args()
+    if opts.idperl is None:
+        opts.idperl = os.path.join(os.getcwd(), "ID_Params.perl")
+    if opts.MA is None:
+        opts.MA = 0.5
+    if opts.MB is None:
+        opts.MB = 0.5
+    opts.t = "bbh"
+
+    # call SpEC script with horizons.h5 and params.perl as input and -t bbh
+    # Import functions from SpEC until we have ported them over. These functions
+    # call old Fortran code (LSODA) through scipy.integrate.odeint, which leads
+    # to lots of noise in stdout. When porting these functions, we should
+    # modernize them to use scipy.integrate.solve_ivp.
+    if "SPEC_ROOT" in os.environ:
+        sys.path.append(os.path.join(os.environ["SPEC_ROOT"], "Support/bin/"))
+
+    try:
+        from OmegaDotEccRemoval import main
+    except ImportError:
+        raise ImportError(
+            "Importing from SpEC failed. Make sure you have pointed "
+            "'-D SPEC_ROOT' to a SpEC installation when configuring the build "
+            "with CMake or set SPEC_ROOT as an environment variable now."
+        )
+
+    editHorizons(opts)
+    main(opts)