Coordinates and TwoPunctures from EOB IDs¶
This tutorial has two goals:
illustrate the main utilities in
PyART.analytic.CoordsChange;show how to start from EOB initial data, map them to ADM with
eob_ID_to_ADM, and assemble a TwoPunctures-ready initial-data dictionary.
The low-level coordinate transforms are useful for inspection and sanity checks. For the actual EOB-to-ADM handoff, the main function to use is eob_ID_to_ADM.
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from pathlib import Path
from tempfile import TemporaryDirectory
import numpy as np
import matplotlib.pyplot as plt
from pprint import pprint
from PyART.analytic import CoordsChange
import PyART.analytic.coordschange as cc
from PyART.numerical.twopuncts_ID import TwoPunctID
try:
import PyART.models.teob as teob
HAVE_TEOB = True
TEOB_IMPORT_ERROR = None
except Exception as exc:
HAVE_TEOB = False
TEOB_IMPORT_ERROR = exc
np.set_printoptions(precision=10, suppress=True)
print(f"TEOB available = {HAVE_TEOB}")
if TEOB_IMPORT_ERROR is not None:
print(f"TEOB import error: {TEOB_IMPORT_ERROR}")
print("To run this tutorial, install TEOBResumS and ensure it is in your PYTHONPATH.")
WARNING: TEOBResumS not installed.
TEOB available = True
1. Low-level coordinate changes¶
CoordsChange exposes polar/cartesian transforms and the 2PN EOB/ADM maps. The small example below performs both round trips and reports the residuals.
r0, phi0, pr0, pphi0 = 12.0, 0.3, -0.02, 4.6
x, y, px, py = CoordsChange.Polar2Cartesian(r0, phi0, pr0, pphi0)
r1, phi1, pr1, pphi1 = CoordsChange.Cartesian2Polar(x, y, px, py)
q = 1.2
nu = q / (1.0 + q) ** 2
qe = np.array([x, y])
pe = np.array([px, py])
qa, pa = CoordsChange.Eob2Adm(qe, pe, nu, PN_order=2)
qe_back, pe_back = CoordsChange.Adm2Eob(qa, pa, nu, PN_order=2)
print("polar -> cartesian")
print(np.array([x, y, px, py]))
print()
print("cartesian -> polar residuals")
print(np.array([r1 - r0, phi1 - phi0, pr1 - pr0, pphi1 - pphi0]))
print()
print("EOB -> ADM")
print("q_ADM =", qa)
print("p_ADM =", pa)
print()
print("round-trip residuals")
print("max |qe_back - qe| =", np.max(np.abs(qe_back - qe)))
print("max |pe_back - pe| =", np.max(np.abs(pe_back - pe)))
polar -> cartesian
[11.4640378695 3.5462424799 -0.1323894757 0.3603019167]
cartesian -> polar residuals
[ 0. -0. -0. 0.]
EOB -> ADM
q_ADM = [10.5861678003 3.2569829327]
p_ADM = [-0.1413954526 0.3907243109]
round-trip residuals
max |qe_back - qe| = 0.00652476010515457
max |pe_back - pe| = 0.000497622413047405
2. Reusable EOB and TwoPunctures helpers¶
The functions do the following:
ADM_from_EOB(...)builds a TEOB waveform, plots the \((2,2)\) amplitude, computes \(t_{\rm mrg} = -u_0\), and converts the initial state to ADM witheob_ID_to_ADM.TP_from_ADM(...)maps that ADM dictionary into the TwoPunctures-style initial-data dictionary.TP_diagnostics_from_ADM(...)keeps a few geometric diagnostics explicit, including the center offset.TP_parfile_from_ADM(...)writes a TwoPunctures parfile from the sameadm_iddata when the configuration is compatible with the currentTwoPunctIDhelper.
def ADM_from_EOB(
M=1.0,
q=1.0,
chi1=(0.0, 0.0, 0.0),
chi2=(0.0, 0.0, 0.0),
f0=0.001,
e0=0.1,
plot_waveform=True,
verbose=True,
):
"""Initial data based on EOB prediction, mapped to ADM."""
if not HAVE_TEOB:
raise RuntimeError(
"TEOBResumS is not available in this environment, so ADM_from_EOB cannot run."
)
pars = teob.CreateDict(
M=M,
q=q,
chi1x=chi1[0],
chi1y=chi1[1],
chi1z=chi1[2],
chi2x=chi2[0],
chi2y=chi2[1],
chi2z=chi2[2],
f0=f0,
ecc=e0,
)
eob_wave = teob.Waveform_EOB(pars=pars)
if plot_waveform:
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(eob_wave.u, eob_wave.hlm[(2, 2)]["A"])
ax.set_xlabel("t / M")
ax.set_ylabel(r"$|h_{22}|$")
ax.set_title("TEOB waveform amplitude")
ax.grid(True, alpha=0.3)
plt.show()
t_mrg = -eob_wave.u[0]
adm_redundant_id = cc.eob_ID_to_ADM(eob_wave, verbose=verbose)
adm_id = {
"px": adm_redundant_id["px"],
"py": adm_redundant_id["py"],
"D": adm_redundant_id["D"],
"M": M,
"q": q,
"chi1x": chi1[0],
"chi1y": chi1[1],
"chi1z": chi1[2],
"chi2x": chi2[0],
"chi2y": chi2[1],
"chi2z": chi2[2],
"x1": adm_redundant_id["x1"],
"x2": adm_redundant_id["x2"],
"x_offset": adm_redundant_id["x_offset"],
}
return adm_id, t_mrg
def TP_from_ADM(adm_id):
"""Build a TwoPunctures-style initial-data dictionary from ADM IDs."""
D = adm_id["D"]
M = adm_id["M"]
q = adm_id["q"]
assert q >= 1.0
mp = M * q / (1 + q)
mm = M / (1 + q)
xp = D * mm
xm = -D * mp
Spx = adm_id["chi1x"] * mp**2
Spy = adm_id["chi1y"] * mp**2
Spz = adm_id["chi1z"] * mp**2
Smx = adm_id["chi2x"] * mm**2
Smy = adm_id["chi2y"] * mm**2
Smz = adm_id["chi2z"] * mm**2
Ppx = adm_id["px"]
Ppy = adm_id["py"]
Ppz = 0.0
Pmx = -Ppx
Pmy = -Ppy
Pmz = -Ppz
initial_data = {
"bh_0_m": mp,
"bh_0_madm": mp,
"bh_0_x": xp,
"bh_0_px": Ppx,
"bh_0_py": Ppy,
"bh_0_sx": Spx,
"bh_0_sy": Spy,
"bh_0_sz": Spz,
"bh_1_m": mm,
"bh_1_madm": mm,
"bh_1_x": xm,
"bh_1_px": Pmx,
"bh_1_py": Pmy,
"bh_1_sx": Smx,
"bh_1_sy": Smy,
"bh_1_sz": Smz,
"give_bare_mass": False,
"bitant": False,
}
return initial_data
def TP_diagnostics_from_ADM(adm_id):
"""Expose a few geometric quantities that are useful for checks."""
D = adm_id["D"]
M = adm_id["M"]
q = adm_id["q"]
mp = M * q / (1 + q)
mm = M / (1 + q)
xp = D * mm
xm = -D * mp
half_D = D * M / 2.0
center_offset = xp - half_D
return {
"xp": xp,
"xm": xm,
"half_D": half_D,
"center_offset": center_offset,
}
def TP_parfile_from_ADM(
adm_id,
outdir,
parfile_name="from_eob_ids.par",
npoints_A=8,
npoints_B=8,
npoints_phi=6,
verbose=False,
):
"""Write a TwoPunctures parfile from ADM IDs using TwoPunctID."""
if abs(adm_id["chi1x"]) > 0 or abs(adm_id["chi1y"]) > 0:
raise ValueError(
"TwoPunctID currently supports only aligned spins in the parfile helper."
)
if abs(adm_id["chi2x"]) > 0 or abs(adm_id["chi2y"]) > 0:
raise ValueError(
"TwoPunctID currently supports only aligned spins in the parfile helper."
)
D = float(adm_id["D"])
px = float(adm_id["px"])
py = float(adm_id["py"])
P = float(np.hypot(px, py))
L = float(D * py)
tp = TwoPunctID(
q=float(adm_id["q"]),
D=D,
L=L,
chi1z=float(adm_id["chi1z"]),
chi2z=float(adm_id["chi2z"]),
npoints_A=npoints_A,
npoints_B=npoints_B,
npoints_phi=npoints_phi,
outdir=outdir,
verbose=verbose,
)
parfile = tp.create_TP_parfile(P, parfile=parfile_name)
return Path(outdir) / parfile
3. Run the workflow on a concrete configuration¶
We now reproduce the same flow as in your driver script: choose intrinsic parameters, call ADM_from_EOB, convert them with TP_from_ADM, and inspect the final TwoPunctures dictionary.
The cell also checks a few consistency relations, including the agreement between the geometric center offset and the x_offset returned by eob_ID_to_ADM.
intrinsic = {
"test": {
"q": 2.0,
"chi1x": 0.0,
"chi1y": 0.0,
"chi1z": 0.0,
"chi2x": 0.0,
"chi2y": 0.0,
"chi2z": 0.0,
"f0": 0.01006,
"ecc": 0.1,
}
}
ID = "test"
adm_id, tmax = ADM_from_EOB(
M=1.0,
q=intrinsic[ID]["q"],
chi1=[
intrinsic[ID]["chi1x"],
intrinsic[ID]["chi1y"],
intrinsic[ID]["chi1z"],
],
chi2=[
intrinsic[ID]["chi2x"],
intrinsic[ID]["chi2y"],
intrinsic[ID]["chi2z"],
],
f0=intrinsic[ID]["f0"],
e0=intrinsic[ID]["ecc"],
plot_waveform=False,
verbose=False,
)
initial_data = TP_from_ADM(adm_id)
tp_diagnostics = TP_diagnostics_from_ADM(adm_id)
assert np.isclose(initial_data["bh_0_x"] - initial_data["bh_1_x"], adm_id["D"])
assert np.isclose(initial_data["bh_0_px"], -initial_data["bh_1_px"])
assert np.isclose(initial_data["bh_0_py"], -initial_data["bh_1_py"])
assert np.isclose(tp_diagnostics["center_offset"], adm_id["x_offset"])
print("ADM IDs")
for key in ("D", "px", "py", "x1", "x2", "x_offset"):
print(f"{key:10s} = {adm_id[key]}")
print(f"{'t_mrg':10s} = {tmax}")
print("\nTwoPunctures initial_data")
for key in (
"bh_0_m",
"bh_0_x",
"bh_0_px",
"bh_0_py",
"bh_0_sz",
"bh_1_m",
"bh_1_x",
"bh_1_px",
"bh_1_py",
"bh_1_sz",
"give_bare_mass",
"bitant",
):
print(f"{key:14s} = {initial_data[key]}")
print("\nDerived diagnostics")
for key, value in tp_diagnostics.items():
print(f"{key:14s} = {value}")
print("\nFull initial_data dictionary")
pprint(initial_data)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 16
12 }
13 }
14
15 ID = "test"
---> 16 adm_id, tmax = ADM_from_EOB(
17 M=1.0,
18 q=intrinsic[ID]["q"],
19 chi1=[
Cell In[3], line 29, in ADM_from_EOB(M, q, chi1, chi2, f0, e0, plot_waveform, verbose)
25 chi2z=chi2[2],
26 f0=f0,
27 ecc=e0,
28 )
---> 29 eob_wave = teob.Waveform_EOB(pars=pars)
30
31 if plot_waveform:
32 fig, ax = plt.subplots(figsize=(8, 4))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/models/teob.py:27, in Waveform_EOB.__init__(self, pars)
25 self.pars = pars
26 self._kind = "EOB"
---> 27 self._run_py()
28 pass
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/models/teob.py:61, in Waveform_EOB._run_py(self)
59 self.domain = "Freq"
60 else:
---> 61 t, hp, hc, hlm, dyn = EOB.EOBRunPy(self.pars)
62 self._u = t
63 hlm_conv = convert_hlm(hlm)
NameError: name 'EOB' is not defined
4. Create a TwoPunctures parfile from the same ADM IDs¶
The initial_data dictionary is useful when you want the puncture data in Python. In addition, the same adm_id can be turned into a TwoPunctures parfile through TwoPunctID.create_TP_parfile.
For the current TwoPunctID helper, this parfile path is restricted to aligned-spin configurations, so the helper raises if the in-plane spin components are nonzero.
with TemporaryDirectory() as tmpdir:
parfile_path = TP_parfile_from_ADM(adm_id, outdir=tmpdir)
selected_keys = {
"par_b",
"par_P_plus1",
"par_P_plus2",
"par_P_minus1",
"par_P_minus2",
"center_offset1",
}
parfile_values = {}
for line in parfile_path.read_text().splitlines():
if "=" not in line or line.lstrip().startswith("#"):
continue
key, value = line.split("=", 1)
key = key.strip()
if key in selected_keys:
parfile_values[key] = float(value)
expected_parfile_values = {
"par_b": adm_id["D"] / 2.0,
"par_P_plus1": adm_id["px"],
"par_P_plus2": adm_id["py"],
"par_P_minus1": -adm_id["px"],
"par_P_minus2": -adm_id["py"],
"center_offset1": adm_id["x_offset"],
}
print(f"parfile written to: {parfile_path}")
print()
for key in sorted(selected_keys):
actual = parfile_values[key]
expected = expected_parfile_values[key]
print(f"{key:14s} actual = {actual: .12e} expected = {expected: .12e}")
assert np.isclose(actual, expected)
parfile written to: /tmp/tmpwgk547lu/from_eob_ids.par
center_offset1 actual = -1.667638744697e+00 expected = -1.667638744697e+00
par_P_minus1 actual = 0.000000000000e+00 expected = -0.000000000000e+00
par_P_minus2 actual = -8.286810492539e-02 expected = -8.286810492539e-02
par_P_plus1 actual = -0.000000000000e+00 expected = 0.000000000000e+00
par_P_plus2 actual = 8.286810492539e-02 expected = 8.286810492539e-02
par_b actual = 5.002916234090e+00 expected = 5.002916234090e+00
Summary¶
In this notebook we have shown how to convert EOB initial data into ADM puncture data, and then into a TwoPunctures initial-data dictionary. We also showed how to write a TwoPunctures parfile from the same ADM data, and verified the key consistency relations between the geometric center offset, the puncture positions, and the parfile entries.
Note: the parfile helper currently follows the capabilities of TwoPunctID, so it assumes that the
punctures start on the x-axis and with no z-component of the momenta (i.e. "center_offset2": 0.0, "center_offset3": 0.0 and "par_P_plus3": 0.0, "par_P_minus3": 0.0)
in the generated parfile.