From ba379b59abccd9379d9e38622b40d5ef7ada08ec Mon Sep 17 00:00:00 2001 From: Max Briel Date: Thu, 26 Feb 2026 23:06:39 +0100 Subject: [PATCH 01/10] unitless ivp change --- posydon/binary_evol/DT/double_CO.py | 245 +++++++++++++++++----------- 1 file changed, 147 insertions(+), 98 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index a73be3299..e7bbf0028 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -1,4 +1,19 @@ -"""Detached evolution for double compact-object binaries.""" +"""Detached evolution for double compact-object binaries. + +The ODEs are non-dimensionalized using the initial orbital separation a₀ +and the characteristic gravitational-wave inspiral timescale + + t₀ = (5/64) c⁵ a₀⁴ / (G³ M m₁ m₂) + +so that the Peters (1964) equations become unit-free with O(1) coefficients: + + dα/dτ = −1 / [α³ (1−e²)^(7/2)] × (1 + 73/24 e² + 37/96 e⁴) + de/dτ = −(19/12) e / [α⁴ (1−e²)^(5/2)] × (1 + 121/304 e²) + +where α = a/a₀ and τ = t/t₀. This avoids the "step size smaller than +machine epsilon" failures that occur when solve_ivp works with physical +units spanning many orders of magnitude. +""" __authors__ = [ @@ -12,9 +27,7 @@ from scipy.integrate import solve_ivp import posydon.utils.constants as constants -from posydon.binary_evol.binarystar import BINARYPROPERTIES from posydon.binary_evol.DT.step_detached import detached_evolution, detached_step -from posydon.binary_evol.singlestar import STARPROPERTIES from posydon.utils.common_functions import ( CO_radius, orbital_period_from_separation, @@ -23,154 +36,190 @@ from posydon.utils.posydonerror import NumericalError -def event(terminal, direction=0): - """Return a helper function to set attributes for solve_ivp events.""" +def _event(terminal, direction=0): + """Return a decorator that sets solve_ivp event attributes.""" def dec(f): f.terminal = terminal f.direction = direction return f return dec + +class _ScaledDenseOutput: + """Wrapper around OdeSolution that maps dimensionless ↔ physical units. + + Parameters + ---------- + sol : OdeSolution + The dense output from solve_ivp in dimensionless variables. + a0_km : float + Initial separation scale factor [km]. + t_scale : float + Characteristic timescale [yr]. + t0_phys : float + Physical time at the start of integration [yr]. + """ + + def __init__(self, sol, a0_km, t_scale, t0_phys): + self.sol = sol + self.a0_km = a0_km + self.t_scale = t_scale + self.t0_phys = t0_phys + + def __call__(self, t_phys): + tau = (np.asarray(t_phys) - self.t0_phys) / self.t_scale + y = self.sol(tau) + y[0] = y[0] * self.a0_km # α → km + return y + + class DoubleCO(detached_step): """Evolve a double compact-object binary due to gravitational radiation. The binary will be evolved until the two compact objects come into contact or until maximum simulation time, based on the quadrupole approximation - of gravitational radiation. - + of gravitational radiation. The integration is performed in dimensionless + variables to guarantee numerical stability. """ def __init__(self, **kwargs): - super().__init__(**kwargs) - - # Reassign detached_evolution to handle - # special DCO evolution + # Use the DCO-specific evolution object (disables tides, winds, etc.) self.evo = double_CO_evolution(**self.evo_kwargs) def __call__(self, binary): - super().__call__(binary) - binary.separation = self.res.y[0][-1] * 100_000 / constants.Rsun + # Convert final separation from km → Rsun + binary.separation = self.res.y[0][-1] * 1e5 / constants.Rsun binary.orbital_period = orbital_period_from_separation( binary.separation, binary.star_1.mass, binary.star_2.mass ) binary.V_sys = binary.V_sys_history[-1] - #if self.res.status == -1: - # set_binary_to_failed(binary) - # raise NumericalError(f"Integration failed for {binary.state} DCO " - # f"({binary.star_1.state}, {binary.star_2.state}): " - # f"{self.res.message}") # contact event triggered if self.res.status == 1: binary.eccentricity = 0.0 binary.state = "contact" binary.event = "CO_contact" - # max time reached + # max time reached (status == 0) elif self.res.status != -1: binary.time = self.max_time binary.eccentricity = self.res.y[1][-1] binary.state = "detached" binary.event = "maxtime" + # ------------------------------------------------------------------ + # Dimensionless ODE integration + # ------------------------------------------------------------------ def solve_ODEs(self, binary, primary, secondary): + """Solve the dimensionless Peters (1964) equations for GW inspiral. + The equations are scaled so that all dependent variables are O(1), + eliminating the "step size < machine epsilon" problem that plagues + the dimensional formulation. A single call to ``solve_ivp`` is + sufficient. + """ self.max_time = binary.properties.max_simulation_time - time_sol = [] - sol = [] - t0 = binary.time - n = 0 - a = binary.separation * constants.Rsun / 100_000 - e = binary.eccentricity - status = -1 - - while status == -1 and n < 6: - try: - res = solve_ivp(self.evo, - events=self.evo.ev_contact, - method="BDF", - t_span=(t0, self.max_time), - y0=[a, e, + t0_phys = binary.time # [yr] + + # --- physical scales ------------------------------------------------ + a0_km = binary.separation * constants.Rsun / 1e5 # [km] + e0 = binary.eccentricity + a0_cgs = a0_km * 1e5 # [cm] + + m1_cgs = primary.mass * constants.msol # [g] + m2_cgs = secondary.mass * constants.msol # [g] + G = constants.standard_cgrav # [cm³ g⁻¹ s⁻²] + c = constants.clight # [cm s⁻¹] + + # characteristic GW inspiral timescale [yr] + t_scale = ((5.0 / 64.0) * c**5 * a0_cgs**4 + / (G**3 * (m1_cgs + m2_cgs) * m1_cgs * m2_cgs) + / constants.secyer) + + # dimensionless contact threshold α_contact = (r1 + r2) / a0 + r1_km = (CO_radius(primary.mass, primary.state) + * constants.Rsun / 1e5) # [km] + r2_km = (CO_radius(secondary.mass, secondary.state) + * constants.Rsun / 1e5) # [km] + alpha_contact = (r1_km + r2_km) / a0_km + + # dimensionless integration limits + tau_max = (self.max_time - t0_phys) / t_scale + + # --- dimensionless RHS (Peters 1964) -------------------------------- + def rhs(tau, y): + alpha, ecc = y[0], y[1] + e2 = ecc * ecc + one_minus_e2 = 1.0 - e2 + + dalpha = (-1.0 + / (alpha**3 * one_minus_e2**3.5) + * (1.0 + (73.0 / 24.0) * e2 + + (37.0 / 96.0) * e2 * e2)) + + decc = (-(19.0 / 12.0) * ecc + / (alpha**4 * one_minus_e2**2.5) + * (1.0 + (121.0 / 304.0) * e2)) + + # ω_sec and ω_pri are constant for DCO (no tides / MB) + return [dalpha, decc, 0.0, 0.0] + + # --- contact event: α hits α_contact from above -------------------- + @_event(True, -1) + def ev_contact(tau, y): + return y[0] - alpha_contact + + # --- integrate ------------------------------------------------------- + try: + res = solve_ivp(rhs, + events=ev_contact, + method="BDF", + t_span=(0.0, tau_max), + y0=[1.0, e0, secondary.omega0, primary.omega0], rtol=1e-10, atol=1e-10, dense_output=True) - - except Exception as e: - set_binary_to_failed(binary) - raise NumericalError(f"SciPy encountered termination edge case while solving GR equations: {e}") - - t0 += res.t[-1] - status = res.status - n += 1 - a = res.y[0][-1] - e = res.y[1][-1] - time_sol.append(res.t) - sol.append(res) + except Exception as exc: + set_binary_to_failed(binary) + raise NumericalError( + "SciPy encountered an error while solving dimensionless " + f"GR equations: {exc}" + ) + + # --- map solution back to physical units ---------------------------- + # Wrap dense output *before* mutating res.t (need original tau range) + res.sol = _ScaledDenseOutput(res.sol, a0_km, t_scale, t0_phys) + + res.y[0] = res.y[0] * a0_km # α → km + res.t = res.t * t_scale + t0_phys # τ → yr + + for i in range(len(res.t_events)): + if len(res.t_events[i]) > 0: + res.t_events[i] = res.t_events[i] * t_scale + t0_phys + for i in range(len(res.y_events)): + if len(res.y_events[i]) > 0: + res.y_events[i][:, 0] *= a0_km return res - class double_CO_evolution(detached_evolution): + """Evolution object for double compact-object binaries. - def __init__(self, **kwargs): + Only gravitational radiation is active; tides, winds, and magnetic + braking are disabled. The actual ODE right-hand side is defined as a + dimensionless local function inside ``DoubleCO.solve_ODEs``; this class + exists so that the ``detached_step`` pipeline (track matching, property + updates) has a valid ``evo`` object. + """ + def __init__(self, **kwargs): super().__init__(**kwargs) - - # For DCO system, only gravitational radiation is considered self.do_magnetic_braking = False self.do_tides = False self.do_wind_loss = False self.do_stellar_evolution_and_spin_from_winds = False self.do_gravitational_radiation = True - - - def set_stars(self, primary, secondary, t0=0.0): - - super().set_stars(primary, secondary, t0) - - # calculate CO radii in kilometers (km rather than Rsol for better scaling in ODE solutions) - self.r1 = CO_radius(self.primary.mass, self.primary.state) * constants.Rsun / 100_000 # [km] - self.r2 = CO_radius(self.secondary.mass, self.secondary.state) * constants.Rsun / 100_000 # [km] - - @event(True, -1) - def ev_contact(self, t, y): - # stop when binary separation = r1+r2 [km] - return y[0] - (self.r1 + self.r2) - - def gravitational_radiation(self): - """Calculate the change in orbital separation and eccentricity due to - gravitational wave radiation from two point masses. Calculated as - according to [1]. - - [1] Peters, P. C. (1964). Gravitational Radiation and the - Motion of Two Point Masses. Physical Review, 136(4B), - B1224–B1232 - """ - g = constants.standard_cgrav - c = constants.clight - - a = self.a - e = self.e - - m1 = self.primary.latest["mass"] * constants.msol - m2 = self.secondary.latest["mass"] * constants.msol - - a = a * 100_000 - - da_gr = ((-64 / 5) - * (g ** 3 * (m1 + m2) * m1 * m2) - / (1 - e ** 2) ** (7 / 2) / a ** 3 / c ** 5 - * (1 + (73 / 24) * e ** 2 + (37 / 96) * e ** 4) - * (constants.secyer / 100_000)) - de_gr = ((-304 / 15) - * e * (g ** 3 * (m1 + m2) * m1 * m2) - / (1 - e ** 2) ** (5 / 2) / a ** 4 / c ** 5 - * (1 + (121 / 304) * e ** 2) - * constants.secyer) - - self.da += da_gr - self.de += de_gr From 07bf2a84905bf44aec5b3f1f1a3213ae923394ec Mon Sep 17 00:00:00 2001 From: Max Briel Date: Fri, 27 Feb 2026 15:21:48 +0100 Subject: [PATCH 02/10] change intergrated equations --- posydon/binary_evol/DT/double_CO.py | 203 ++++++++++++++++++---------- 1 file changed, 128 insertions(+), 75 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index e7bbf0028..1b97958e0 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -1,18 +1,19 @@ """Detached evolution for double compact-object binaries. -The ODEs are non-dimensionalized using the initial orbital separation a₀ -and the characteristic gravitational-wave inspiral timescale +The Peters (1964) equations are solved using the substitution s = -ln(alpha), +where α = a/a₀ is the dimensionless separation. This eliminates the +singular α⁻³ and α⁻⁴ prefactors that cause "step size smaller than +machine epsilon" failures for tight compact-object binaries (especially +BH-BH systems whose contact separation α_contact is extremely small). - t₀ = (5/64) c⁵ a₀⁴ / (G³ M m₁ m₂) +With this substitution, the Peters equations become: -so that the Peters (1964) equations become unit-free with O(1) coefficients: + de/ds = -(19/12) e (1-e^2)(1 + 121/304 e^2) / F(e) + dτ/ds = exp(-4s) (1-e^2)^(7/2) / F(e) - dα/dτ = −1 / [α³ (1−e²)^(7/2)] × (1 + 73/24 e² + 37/96 e⁴) - de/dτ = −(19/12) e / [α⁴ (1−e²)^(5/2)] × (1 + 121/304 e²) - -where α = a/a₀ and τ = t/t₀. This avoids the "step size smaller than -machine epsilon" failures that occur when solve_ivp works with physical -units spanning many orders of magnitude. +where F(e) = 1 + 73/24 e² + 37/96 e⁴, and the characteristic timescale +is t₀ = (5/64) c^5 a0^4 / (G^3 M m_1 m_2). The eccentricity ODE is now +autonomous (independent of s), and the system is non-stiff. """ @@ -25,6 +26,7 @@ import numpy as np from scipy.integrate import solve_ivp +from scipy.interpolate import PchipInterpolator import posydon.utils.constants as constants from posydon.binary_evol.DT.step_detached import detached_evolution, detached_step @@ -45,32 +47,63 @@ def dec(f): return dec -class _ScaledDenseOutput: - """Wrapper around OdeSolution that maps dimensionless ↔ physical units. +class _SLogDenseOutput: + """Dense output for the s = -ln(alpha) formulation of the Peters equations. + + Inverts the monotone τ(s) mapping via a PCHIP interpolant, then + evaluates the native ODE dense output in s-space and converts back + to physical variables [separation_Rsun, eccentricity, ω_sec, ω_pri]. Parameters ---------- sol : OdeSolution - The dense output from solve_ivp in dimensionless variables. - a0_km : float - Initial separation scale factor [km]. + Dense output from solve_ivp in s-space. + a0_Rsun : float + Initial separation [Rsun]. t_scale : float - Characteristic timescale [yr]. + Characteristic GW timescale [yr]. t0_phys : float - Physical time at the start of integration [yr]. + Physical time at integration start [yr]. + s_vals : array + Independent-variable values produced by the solver. + tau_vals : array + Corresponding dimensionless-time values τ(s). """ - def __init__(self, sol, a0_km, t_scale, t0_phys): + def __init__(self, sol, a0_Rsun, t_scale, t0_phys, s_vals, tau_vals): self.sol = sol - self.a0_km = a0_km + self.a0_Rsun = a0_Rsun self.t_scale = t_scale self.t0_phys = t0_phys + # Near contact dtau/ds ≈ exp(−4s) → 0, so consecutive tau values + # can be identical to machine precision. Keep only points where + # tau strictly increases so that PchipInterpolator is happy. + tau_vals = np.asarray(tau_vals, dtype=float) + s_vals = np.asarray(s_vals, dtype=float) + mask = np.concatenate(([True], np.diff(tau_vals) > 0)) + tau_mono = tau_vals[mask] + s_mono = s_vals[mask] + if len(tau_mono) >= 2: + self._tau_to_s = PchipInterpolator(tau_mono, s_mono) + else: + # Fallback: constant (binary barely evolved) + _s0 = s_mono[0] if len(s_mono) else 0.0 + self._tau_to_s = lambda tau: np.full_like( + np.asarray(tau, dtype=float), _s0) def __call__(self, t_phys): + # convert to dimensionless time tau = (np.asarray(t_phys) - self.t0_phys) / self.t_scale - y = self.sol(tau) - y[0] = y[0] * self.a0_km # α → km - return y + # convert to s-space + s = np.asarray(self._tau_to_s(tau)) + raw = self.sol(s) # [ecc, tau, secondary.omega, primary.omega] + result = np.empty_like(raw) + alpha = np.exp(-s) + result[0] = alpha * self.a0_Rsun # separation (Rsun) + result[1] = raw[0] # eccentricity + result[2] = raw[2] # omega_sec + result[3] = raw[3] # omega_pri + return result class DoubleCO(detached_step): @@ -90,94 +123,90 @@ def __init__(self, **kwargs): def __call__(self, binary): super().__call__(binary) - # Convert final separation from km → Rsun - binary.separation = self.res.y[0][-1] * 1e5 / constants.Rsun - binary.orbital_period = orbital_period_from_separation( - binary.separation, binary.star_1.mass, binary.star_2.mass - ) binary.V_sys = binary.V_sys_history[-1] - # contact event triggered if self.res.status == 1: binary.eccentricity = 0.0 binary.state = "contact" binary.event = "CO_contact" - # max time reached (status == 0) + elif self.res.status != -1: binary.time = self.max_time binary.eccentricity = self.res.y[1][-1] binary.state = "detached" binary.event = "maxtime" - # ------------------------------------------------------------------ - # Dimensionless ODE integration - # ------------------------------------------------------------------ def solve_ODEs(self, binary, primary, secondary): - """Solve the dimensionless Peters (1964) equations for GW inspiral. - - The equations are scaled so that all dependent variables are O(1), - eliminating the "step size < machine epsilon" problem that plagues - the dimensional formulation. A single call to ``solve_ivp`` is - sufficient. + """Solve the Peters (1964) GW inspiral. + + We do two transformations to the standard Peters equations: + 1. We remove the dimensionality by using the dimensionless separation \alpha = a/a0 and + dimensionless time \tau = t / t0, where a0 is the initial separation + and t0 is the characteristic GW timescale. + 2. We then substitute s = -ln(\alpha) to eliminate the singularity in the Peters equations, + which removes the alpha^-3 and alpha^-4 prefactors that cause numerical issues for tight binaries. """ self.max_time = binary.properties.max_simulation_time t0_phys = binary.time # [yr] # --- physical scales ------------------------------------------------ - a0_km = binary.separation * constants.Rsun / 1e5 # [km] + a0_Rsun = binary.separation # [Rsun] + a0_km = a0_Rsun * constants.Rsun / 1e5 # [km] e0 = binary.eccentricity a0_cgs = a0_km * 1e5 # [cm] m1_cgs = primary.mass * constants.msol # [g] m2_cgs = secondary.mass * constants.msol # [g] - G = constants.standard_cgrav # [cm³ g⁻¹ s⁻²] - c = constants.clight # [cm s⁻¹] + G = constants.standard_cgrav # [cm^3 g^-1 s^02] + c = constants.clight # [cm s^-1] # characteristic GW inspiral timescale [yr] t_scale = ((5.0 / 64.0) * c**5 * a0_cgs**4 / (G**3 * (m1_cgs + m2_cgs) * m1_cgs * m2_cgs) / constants.secyer) - # dimensionless contact threshold α_contact = (r1 + r2) / a0 + # dimensionless contact threshold \alpha_contact = (r1 + r2) / a0 r1_km = (CO_radius(primary.mass, primary.state) * constants.Rsun / 1e5) # [km] r2_km = (CO_radius(secondary.mass, secondary.state) * constants.Rsun / 1e5) # [km] + + # Unitless separation at contact. alpha_contact = (r1_km + r2_km) / a0_km - # dimensionless integration limits + # dimensionless max-time limit tau_max = (self.max_time - t0_phys) / t_scale - # --- dimensionless RHS (Peters 1964) -------------------------------- - def rhs(tau, y): - alpha, ecc = y[0], y[1] + s_contact = -np.log(alpha_contact) + + def rhs(s, y): + ecc = y[0] e2 = ecc * ecc one_minus_e2 = 1.0 - e2 - dalpha = (-1.0 - / (alpha**3 * one_minus_e2**3.5) - * (1.0 + (73.0 / 24.0) * e2 - + (37.0 / 96.0) * e2 * e2)) + f_e = 1.0 + (73.0 / 24.0) * e2 + (37.0 / 96.0) * e2 * e2 + g_e = 1.0 + (121.0 / 304.0) * e2 - decc = (-(19.0 / 12.0) * ecc - / (alpha**4 * one_minus_e2**2.5) - * (1.0 + (121.0 / 304.0) * e2)) + decc_ds = -(19.0 / 12.0) * ecc * one_minus_e2 * g_e / f_e + dtau_ds = np.exp(-4.0 * s) * one_minus_e2**3.5 / f_e - # ω_sec and ω_pri are constant for DCO (no tides / MB) - return [dalpha, decc, 0.0, 0.0] + # primary.omega0 and secondary.omega0 are constant for DCO (no tides / MB) + return [decc_ds, dtau_ds, 0.0, 0.0] - # --- contact event: α hits α_contact from above -------------------- - @_event(True, -1) - def ev_contact(tau, y): - return y[0] - alpha_contact + # --- max-time event: tau reaches tau_max before contact --------------- + @_event(True, +1) + def ev_maxtime(s, y): + return y[1] - tau_max # --- integrate ------------------------------------------------------- + # State vector: [e, tau, secondary.omega, primary.omega] + # Independent variable: s = −ln(a/a0), from 0 to s_contact try: res = solve_ivp(rhs, - events=ev_contact, - method="BDF", - t_span=(0.0, tau_max), - y0=[1.0, e0, + events=ev_maxtime, + method="RK45", + t_span=(0.0, s_contact), + y0=[e0, 0.0, secondary.omega0, primary.omega0], rtol=1e-10, atol=1e-10, @@ -185,23 +214,47 @@ def ev_contact(tau, y): except Exception as exc: set_binary_to_failed(binary) raise NumericalError( - "SciPy encountered an error while solving dimensionless " - f"GR equations: {exc}" + "SciPy encountered an error while solving " + f"GR equations (s-formulation): {exc}" ) # --- map solution back to physical units ---------------------------- - # Wrap dense output *before* mutating res.t (need original tau range) - res.sol = _ScaledDenseOutput(res.sol, a0_km, t_scale, t0_phys) - - res.y[0] = res.y[0] * a0_km # α → km - res.t = res.t * t_scale + t0_phys # τ → yr - + s_vals = res.t.copy() + alpha_vals = np.exp(-s_vals) + ecc_vals = res.y[0].copy() + tau_vals = res.y[1].copy() + + # Wrap dense output to convert from s-space back to physical variables + res.sol = _SLogDenseOutput(res.sol, a0_Rsun, t_scale, t0_phys, + s_vals, tau_vals) + + # Map solver status to the convention expected by __call__: + # status 0 -> reached s_contact -> contact -> report as 1 + # status 1 -> maxtime event -> no merge -> report as 0 + # status −1 -> solver failure -> keep as −1 + if res.status == 0: + res.status = 1 # contact + elif res.status == 1: + res.status = 0 # maxtime + + # Rearrange state vector from [ecc, tau, secondary.omega, primary.omega] in s-space + # to [sep_Rsun, ecc, secondary.omega, primary.omega] + res.y[0] = alpha_vals * a0_Rsun # sep (Rsun) + res.y[1] = ecc_vals # ecc + # res.y[2], res.y[3] unchanged (omegas are constant for DCO) + + res.t = tau_vals * t_scale + t0_phys # tau -> yr + + # Convert events from s-space to physical units for i in range(len(res.t_events)): if len(res.t_events[i]) > 0: - res.t_events[i] = res.t_events[i] * t_scale + t0_phys - for i in range(len(res.y_events)): - if len(res.y_events[i]) > 0: - res.y_events[i][:, 0] *= a0_km + s_ev = res.t_events[i].copy() + alpha_ev = np.exp(-s_ev) + ecc_ev = res.y_events[i][:, 0].copy() + tau_ev = res.y_events[i][:, 1].copy() + res.y_events[i][:, 0] = alpha_ev * a0_Rsun # sep (Rsun) + res.y_events[i][:, 1] = ecc_ev # ecc + res.t_events[i] = tau_ev * t_scale + t0_phys # yr return res From 68c114caf515a8c723762d91208a06cce14e2f54 Mon Sep 17 00:00:00 2001 From: Max Briel Date: Mon, 2 Mar 2026 14:38:03 +0100 Subject: [PATCH 03/10] cleanup --- posydon/binary_evol/DT/double_CO.py | 45 +++++++++++++++++------------ 1 file changed, 27 insertions(+), 18 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index 1b97958e0..f816266b4 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -75,7 +75,7 @@ def __init__(self, sol, a0_Rsun, t_scale, t0_phys, s_vals, tau_vals): self.a0_Rsun = a0_Rsun self.t_scale = t_scale self.t0_phys = t0_phys - # Near contact dtau/ds ≈ exp(−4s) → 0, so consecutive tau values + # Near contact dtau/ds ≈ exp(−4s) -> 0, so consecutive tau values # can be identical to machine precision. Keep only points where # tau strictly increases so that PchipInterpolator is happy. tau_vals = np.asarray(tau_vals, dtype=float) @@ -179,20 +179,6 @@ def solve_ODEs(self, binary, primary, secondary): s_contact = -np.log(alpha_contact) - def rhs(s, y): - ecc = y[0] - e2 = ecc * ecc - one_minus_e2 = 1.0 - e2 - - f_e = 1.0 + (73.0 / 24.0) * e2 + (37.0 / 96.0) * e2 * e2 - g_e = 1.0 + (121.0 / 304.0) * e2 - - decc_ds = -(19.0 / 12.0) * ecc * one_minus_e2 * g_e / f_e - dtau_ds = np.exp(-4.0 * s) * one_minus_e2**3.5 / f_e - - # primary.omega0 and secondary.omega0 are constant for DCO (no tides / MB) - return [decc_ds, dtau_ds, 0.0, 0.0] - # --- max-time event: tau reaches tau_max before contact --------------- @_event(True, +1) def ev_maxtime(s, y): @@ -237,11 +223,11 @@ def ev_maxtime(s, y): elif res.status == 1: res.status = 0 # maxtime - # Rearrange state vector from [ecc, tau, secondary.omega, primary.omega] in s-space - # to [sep_Rsun, ecc, secondary.omega, primary.omega] + # Rearrange state vector in s-space + # from [ecc, tau, secondary.omega, primary.omega] + # to [sep_Rsun, ecc, secondary.omega, primary.omega] res.y[0] = alpha_vals * a0_Rsun # sep (Rsun) res.y[1] = ecc_vals # ecc - # res.y[2], res.y[3] unchanged (omegas are constant for DCO) res.t = tau_vals * t_scale + t0_phys # tau -> yr @@ -276,3 +262,26 @@ def __init__(self, **kwargs): self.do_wind_loss = False self.do_stellar_evolution_and_spin_from_winds = False self.do_gravitational_radiation = True + + def rhs(s, y): + """Right-hand side of the ODE system in s-space. + + Parameters + ---------- + s : float + Independent variable, s = -ln(a/a0). + y : array_like + State vector at s, [ecc, tau, secondary.omega, primary.omega]. + """ + ecc = y[0] + e2 = ecc * ecc + one_minus_e2 = 1.0 - e2 + + f_e = 1.0 + (73.0 / 24.0) * e2 + (37.0 / 96.0) * e2 * e2 + g_e = 1.0 + (121.0 / 304.0) * e2 + + decc_ds = -(19.0 / 12.0) * ecc * one_minus_e2 * g_e / f_e + dtau_ds = np.exp(-4.0 * s) * one_minus_e2**3.5 / f_e + + # primary.omega0 and secondary.omega0 are constant for DCO (no tides / MB) + return [decc_ds, dtau_ds, 0.0, 0.0] From 497a4c89f58c2c32b92b06274910f59fdd8706bc Mon Sep 17 00:00:00 2001 From: Max Briel Date: Mon, 2 Mar 2026 15:00:58 +0100 Subject: [PATCH 04/10] remove description --- posydon/binary_evol/DT/double_CO.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index 68e68f027..5090cc6c6 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -1,19 +1,4 @@ """Detached evolution for double compact-object binaries. - -The Peters (1964) equations are solved using the substitution s = -ln(alpha), -where α = a/a₀ is the dimensionless separation. This eliminates the -singular α⁻³ and α⁻⁴ prefactors that cause "step size smaller than -machine epsilon" failures for tight compact-object binaries (especially -BH-BH systems whose contact separation α_contact is extremely small). - -With this substitution, the Peters equations become: - - de/ds = -(19/12) e (1-e^2)(1 + 121/304 e^2) / F(e) - dτ/ds = exp(-4s) (1-e^2)^(7/2) / F(e) - -where F(e) = 1 + 73/24 e² + 37/96 e⁴, and the characteristic timescale -is t₀ = (5/64) c^5 a0^4 / (G^3 M m_1 m_2). The eccentricity ODE is now -autonomous (independent of s), and the system is non-stiff. """ From b419d227ff0ec04077dbdf90b31d09114eebf476 Mon Sep 17 00:00:00 2001 From: Max Briel Date: Mon, 2 Mar 2026 15:01:46 +0100 Subject: [PATCH 05/10] =?UTF-8?q?remove=20=CF=84=20for=20tau?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- posydon/binary_evol/DT/double_CO.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index 5090cc6c6..e71408c33 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -35,7 +35,7 @@ def dec(f): class _SLogDenseOutput: """Dense output for the s = -ln(alpha) formulation of the Peters equations. - Inverts the monotone τ(s) mapping via a PCHIP interpolant, then + Inverts the monotone tau(s) mapping via a PCHIP interpolant, then evaluates the native ODE dense output in s-space and converts back to physical variables [separation_Rsun, eccentricity, ω_sec, ω_pri]. @@ -52,7 +52,7 @@ class _SLogDenseOutput: s_vals : array Independent-variable values produced by the solver. tau_vals : array - Corresponding dimensionless-time values τ(s). + Corresponding dimensionless-time values tau(s). """ def __init__(self, sol, a0_Rsun, t_scale, t0_phys, s_vals, tau_vals): From d22d12ca90366a260d0e8cadddd0c7c7e0fd3ed8 Mon Sep 17 00:00:00 2001 From: Max Briel Date: Mon, 2 Mar 2026 15:04:08 +0100 Subject: [PATCH 06/10] fix rhs call after moved to evolution class --- posydon/binary_evol/DT/double_CO.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index e71408c33..05993b2eb 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -177,7 +177,7 @@ def ev_maxtime(s, y): # State vector: [e, tau, secondary.omega, primary.omega] # Independent variable: s = −ln(a/a0), from 0 to s_contact try: - res = solve_ivp(rhs, + res = solve_ivp(self.evo.rhs, events=ev_maxtime, method="RK45", t_span=(0.0, s_contact), From df207d03daee19a1574371a78919591c5125dd4a Mon Sep 17 00:00:00 2001 From: Max <14039563+maxbriel@users.noreply.github.com> Date: Thu, 12 Mar 2026 15:25:14 +0100 Subject: [PATCH 07/10] Fix method signature for rhs to include self --- posydon/binary_evol/DT/double_CO.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index 05993b2eb..d75cd6357 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -252,7 +252,7 @@ def __init__(self, **kwargs): self.do_stellar_evolution_and_spin_from_winds = False self.do_gravitational_radiation = True - def rhs(s, y): + def rhs(self, s, y): """Right-hand side of the ODE system in s-space. Parameters From 799d180991b597c327e191aee558e4ec7b2792ba Mon Sep 17 00:00:00 2001 From: Max Briel Date: Thu, 12 Mar 2026 16:13:59 +0100 Subject: [PATCH 08/10] remove call to time (which is now set in detached.call --- posydon/binary_evol/DT/double_CO.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index d75cd6357..78c867c79 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -111,10 +111,6 @@ def __call__(self, binary): binary.V_sys = binary.V_sys_history[-1] if self.res.status == 1: - binary.time = self.res.real_time[-1] - if len(self.res.real_time) > 1: - for k, real_time in enumerate(self.res.real_time[::-1]): - binary.time_history[-(k+1)] = real_time binary.eccentricity = 0.0 binary.state = "contact" binary.event = "CO_contact" From 07bd72fd8cd9daf070e62f2065baaf99eea34235 Mon Sep 17 00:00:00 2001 From: Max Briel Date: Fri, 20 Mar 2026 10:35:59 +0100 Subject: [PATCH 09/10] add eccentricity in ln-space --- posydon/binary_evol/DT/double_CO.py | 129 +++++++++++++++++----------- 1 file changed, 78 insertions(+), 51 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index 78c867c79..708151ec3 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -11,7 +11,7 @@ import numpy as np from scipy.integrate import solve_ivp -from scipy.interpolate import PchipInterpolator +from scipy.optimize import brentq import posydon.utils.constants as constants from posydon.binary_evol.DT.step_detached import detached_evolution, detached_step @@ -32,12 +32,12 @@ def dec(f): return dec -class _SLogDenseOutput: +class _SBrentqDenseOutput: """Dense output for the s = -ln(alpha) formulation of the Peters equations. - Inverts the monotone tau(s) mapping via a PCHIP interpolant, then - evaluates the native ODE dense output in s-space and converts back - to physical variables [separation_Rsun, eccentricity, ω_sec, ω_pri]. + Inverts the monotone tau(s) mapping via ``brentq`` root-finding on the + native ODE dense output, then converts back to physical variables + [separation_Rsun, eccentricity, ω_sec, ω_pri]. Parameters ---------- @@ -49,45 +49,41 @@ class _SLogDenseOutput: Characteristic GW timescale [yr]. t0_phys : float Physical time at integration start [yr]. - s_vals : array - Independent-variable values produced by the solver. - tau_vals : array - Corresponding dimensionless-time values tau(s). + s_lo : float + Lower bound of s integration range. + s_hi : float + Upper bound of s integration range. """ - def __init__(self, sol, a0_Rsun, t_scale, t0_phys, s_vals, tau_vals): + def __init__(self, sol, a0_Rsun, t_scale, t0_phys, s_lo, s_hi): self.sol = sol self.a0_Rsun = a0_Rsun self.t_scale = t_scale self.t0_phys = t0_phys - # Near contact dtau/ds ≈ exp(−4s) -> 0, so consecutive tau values - # can be identical to machine precision. Keep only points where - # tau strictly increases so that PchipInterpolator is happy. - tau_vals = np.asarray(tau_vals, dtype=float) - s_vals = np.asarray(s_vals, dtype=float) - mask = np.concatenate(([True], np.diff(tau_vals) > 0)) - tau_mono = tau_vals[mask] - s_mono = s_vals[mask] - if len(tau_mono) >= 2: - self._tau_to_s = PchipInterpolator(tau_mono, s_mono) - else: - # Fallback: constant (binary barely evolved) - _s0 = s_mono[0] if len(s_mono) else 0.0 - self._tau_to_s = lambda tau: np.full_like( - np.asarray(tau, dtype=float), _s0) + self.s_lo = float(s_lo) + self.s_hi = float(s_hi) def __call__(self, t_phys): - # convert to dimensionless time - tau = (np.asarray(t_phys) - self.t0_phys) / self.t_scale - # convert to s-space - s = np.asarray(self._tau_to_s(tau)) - raw = self.sol(s) # [ecc, tau, secondary.omega, primary.omega] - result = np.empty_like(raw) - alpha = np.exp(-s) - result[0] = alpha * self.a0_Rsun # separation (Rsun) - result[1] = raw[0] # eccentricity - result[2] = raw[2] # omega_sec - result[3] = raw[3] # omega_pri + t_phys = np.atleast_1d(np.asarray(t_phys, dtype=float)) + tau_target = (t_phys - self.t0_phys) / self.t_scale + + n = len(tau_target) + result = np.empty((4, n)) + for i in range(n): + s_star = brentq( + lambda s: self.sol(s)[1] - tau_target[i], + self.s_lo, self.s_hi, + xtol=1e-14, rtol=1e-14, + ) + raw = self.sol(s_star) # [l, tau, omega_sec, omega_pri] + alpha = np.exp(-s_star) + result[0, i] = alpha * self.a0_Rsun # separation (Rsun) + result[1, i] = np.exp(raw[0]) # eccentricity = exp(l) + result[2, i] = raw[2] # omega_sec + result[3, i] = raw[3] # omega_pri + + if result.shape[1] == 1: + return result[:, 0] return result @@ -170,14 +166,17 @@ def ev_maxtime(s, y): return y[1] - tau_max # --- integrate ------------------------------------------------------- - # State vector: [e, tau, secondary.omega, primary.omega] + # State vector: [l, tau, secondary.omega, primary.omega] + # where l = ln(e) and tau = t/t0 # Independent variable: s = −ln(a/a0), from 0 to s_contact + l0 = (np.log(e0) if e0 > 0 + else np.log(np.finfo(float).tiny)) try: res = solve_ivp(self.evo.rhs, events=ev_maxtime, method="RK45", t_span=(0.0, s_contact), - y0=[e0, 0.0, + y0=[l0, 0.0, secondary.omega0, primary.omega0], rtol=1e-10, atol=1e-10, @@ -192,12 +191,13 @@ def ev_maxtime(s, y): # --- map solution back to physical units ---------------------------- s_vals = res.t.copy() alpha_vals = np.exp(-s_vals) - ecc_vals = res.y[0].copy() + l_vals = res.y[0].copy() + ecc_vals = np.exp(l_vals) tau_vals = res.y[1].copy() - # Wrap dense output to convert from s-space back to physical variables - res.sol = _SLogDenseOutput(res.sol, a0_Rsun, t_scale, t0_phys, - s_vals, tau_vals) + # Wrap dense output to convert from ln-space back to physical variables + res.sol = _SBrentqDenseOutput(res.sol, a0_Rsun, t_scale, t0_phys, + s_vals[0], s_vals[-1]) # Map solver status to the convention expected by __call__: # status 0 -> reached s_contact -> contact -> report as 1 @@ -208,8 +208,8 @@ def ev_maxtime(s, y): elif res.status == 1: res.status = 0 # maxtime - # Rearrange state vector in s-space - # from [ecc, tau, secondary.omega, primary.omega] + # Rearrange state vector in ln-space + # from [l, tau, secondary.omega, primary.omega] # to [sep_Rsun, ecc, secondary.omega, primary.omega] res.y[0] = alpha_vals * a0_Rsun # sep (Rsun) res.y[1] = ecc_vals # ecc @@ -221,7 +221,8 @@ def ev_maxtime(s, y): if len(res.t_events[i]) > 0: s_ev = res.t_events[i].copy() alpha_ev = np.exp(-s_ev) - ecc_ev = res.y_events[i][:, 0].copy() + l_ev = res.y_events[i][:, 0].copy() + ecc_ev = np.exp(l_ev) tau_ev = res.y_events[i][:, 1].copy() res.y_events[i][:, 0] = alpha_ev * a0_Rsun # sep (Rsun) res.y_events[i][:, 1] = ecc_ev # ecc @@ -248,25 +249,51 @@ def __init__(self, **kwargs): self.do_stellar_evolution_and_spin_from_winds = False self.do_gravitational_radiation = True + def __call__(self, s, y): + """Differential equation describing the orbital evolution of a double + compact object binary. + + + Parameters + ---------- + s : float + Independent variable, s = -ln(a/a0). + y : array_like + State vector at s, [l, tau, secondary.omega, primary.omega], + where l = ln(e) and tau = t / t0. + """ + + if self.do_gravitational_radiation: + return self.rhs(s, y) + + + + def rhs(self, s, y): - """Right-hand side of the ODE system in s-space. + """Right-hand side of the ODE system in ln-space. Parameters ---------- s : float Independent variable, s = -ln(a/a0). y : array_like - State vector at s, [ecc, tau, secondary.omega, primary.omega]. + State vector at s, [l, tau, secondary.omega, primary.omega], + where l = ln(e). """ - ecc = y[0] - e2 = ecc * ecc + l = y[0] + e = np.exp(l) + e2 = e * e + + if e2 >= 1.0: + return [0.0, 0.0, 0.0, 0.0] + one_minus_e2 = 1.0 - e2 f_e = 1.0 + (73.0 / 24.0) * e2 + (37.0 / 96.0) * e2 * e2 g_e = 1.0 + (121.0 / 304.0) * e2 - decc_ds = -(19.0 / 12.0) * ecc * one_minus_e2 * g_e / f_e + dl_ds = -(19.0 / 12.0) * one_minus_e2 * g_e / f_e dtau_ds = np.exp(-4.0 * s) * one_minus_e2**3.5 / f_e # primary.omega0 and secondary.omega0 are constant for DCO (no tides / MB) - return [decc_ds, dtau_ds, 0.0, 0.0] + return [dl_ds, dtau_ds, 0.0, 0.0] From 7a89d6df5d44aa95b6688882f65e845bffa66915 Mon Sep 17 00:00:00 2001 From: Max Briel Date: Fri, 20 Mar 2026 10:40:39 +0100 Subject: [PATCH 10/10] make consistent with notes --- posydon/binary_evol/DT/double_CO.py | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/posydon/binary_evol/DT/double_CO.py b/posydon/binary_evol/DT/double_CO.py index 708151ec3..58032604b 100644 --- a/posydon/binary_evol/DT/double_CO.py +++ b/posydon/binary_evol/DT/double_CO.py @@ -141,7 +141,8 @@ def solve_ODEs(self, binary, primary, secondary): G = constants.standard_cgrav # [cm^3 g^-1 s^02] c = constants.clight # [cm s^-1] - # characteristic GW inspiral timescale [yr] + # Characteristic circular inspiral timescale t0 = a0^4 / (beta) [yr], + # where beta = (64/5) * G^3 * m1 * m2 * M / c^5. t_scale = ((5.0 / 64.0) * c**5 * a0_cgs**4 / (G**3 * (m1_cgs + m2_cgs) * m1_cgs * m2_cgs) / constants.secyer) @@ -267,7 +268,22 @@ def __call__(self, s, y): return self.rhs(s, y) + def _g(e2): + """Peters (1964) f(e): denominator function for the da/dt equation. + Defined as f(e) = 1 + (73/24)e^2 + (37/96)e^4 in Peters (1964), Eq. 5.11. + Named _g here to match the GW-integration convention. + """ + return 1.0 + (73.0 / 24.0) * e2 + (37.0 / 96.0) * e2 * e2 + + + def _f(e2): + """Peters (1964) g(e): numerator function for the de/dt equation. + + Defined as g(e) = 1 + (121/304)e^2 in Peters (1964), Eq. 5.13. + Named _f here to match the GW-integration convention. + """ + return 1.0 + (121.0 / 304.0) * e2 def rhs(self, s, y): """Right-hand side of the ODE system in ln-space. @@ -288,12 +304,11 @@ def rhs(self, s, y): return [0.0, 0.0, 0.0, 0.0] one_minus_e2 = 1.0 - e2 + G = self._g(e2) + F = self._f(e2) - f_e = 1.0 + (73.0 / 24.0) * e2 + (37.0 / 96.0) * e2 * e2 - g_e = 1.0 + (121.0 / 304.0) * e2 - - dl_ds = -(19.0 / 12.0) * one_minus_e2 * g_e / f_e - dtau_ds = np.exp(-4.0 * s) * one_minus_e2**3.5 / f_e + dl_ds = -(19.0 / 12.0) * one_minus_e2 * F / G + dtau_ds = np.exp(-4.0 * s) * one_minus_e2**3.5 / G # primary.omega0 and secondary.omega0 are constant for DCO (no tides / MB) return [dl_ds, dtau_ds, 0.0, 0.0]