Debugging POSYDON Binary Population Synthesis πο
If you havenβt done so yet, export the path POSYDON environment variables. Set these parameters in your .bash_profile or .zshrc if you use POSYDON regularly.
[ ]:
%env PATH_TO_POSYDON=/Users/simone/Google Drive/github/POSYDON-public/
%env PATH_TO_POSYDON_DATA=/Volumes/T7/
To debug the evolution of a binary systems you should load the POSYDON BinaryPopulation object in memory and then call the restore and evolve. This will run the evolution of the binary system and display any error traceback. You can debug the code by setting breakpoints in the POSYDON code until you devug the specific binary evolution.
The trick is to do this is not dump the binary population to disk but to keep them in memory. This is done by setting the optimize_ram=False.
[ ]:
from posydon.popsyn.binarypopulation import BinaryPopulation
from posydon.binary_evol.simulationproperties import SimulationProperties
from posydon.binary_evol.flow_chart import flow_chart
from posydon.binary_evol.MESA.step_mesa import CO_HeMS_step, MS_MS_step, CO_HMS_RLO_step, CO_HeMS_RLO_step
from posydon.binary_evol.DT.step_detached import detached_step
from posydon.binary_evol.DT.step_disrupted import DisruptedStep
from posydon.binary_evol.DT.step_merged import MergedStep
from posydon.binary_evol.CE.step_CEE import StepCEE
from posydon.binary_evol.SN.step_SN import StepSN
from posydon.binary_evol.DT.double_CO import DoubleCO
from posydon.binary_evol.step_end import step_end
from posydon.binary_evol.simulationproperties import StepNamesHooks
# STEP CUSTOMISATION
MESA_STEP = dict(
interpolation_path = None, # found by default
interpolation_filename = None, # found by default
interpolation_method = 'nearest_neighbour', # 'nearest_neighbour' 'linear3c_kNN' '1NN_1NN'
save_initial_conditions = True, # only for interpolation_method='nearest_neighbour'
track_interpolation = False, # True False
stop_method = 'stop_at_max_time', # 'stop_at_end' 'stop_at_max_time' 'stop_at_condition'
stop_star = 'star_1', # only for stop_method='stop_at_condition' 'star_1' 'star_2'
stop_var_name = None, # only for stop_method='stop_at_condition' str
stop_value = None, # only for stop_method='stop_at_condition' float
stop_interpolate = True, # True False
verbose = False, # True False
)
DETACHED_STEP = dict(
matching_method = 'minimize', #'minimize' 'root'
do_wind_loss = True, # True False
do_tides = True, # True False
do_gravitational_radiation = True, # True False
do_magnetic_braking = True, # True False
do_stellar_evolution_and_spin_from_winds = True, # True False
RLO_orbit_at_orbit_with_same_am = False, # True False
verbose = False, # True False
)
DISRUPTED_STEP = dict(
grid_name_Hrich=None,
grid_name_strippedHe=None,
metallicity=None,
dt=None,
n_o_steps_history=None,
matching_method="minimize",
initial_mass=None,
rootm=None,
verbose=False,
do_wind_loss=True,
do_tides=True,
do_gravitational_radiation=True,
do_magnetic_braking=True,
magnetic_braking_mode="RVJ83",
do_stellar_evolution_and_spin_from_winds=True,
RLO_orbit_at_orbit_with_same_am=False,
list_for_matching_HMS=None,
list_for_matching_postMS=None,
list_for_matching_HeStar=None
)
CE_STEP = dict(
prescription='alpha-lambda', # 'alpha-lambda'
common_envelope_efficiency=1.0, # float in (0, inf)
common_envelope_option_for_lambda='lambda_from_grid_final_values', # (1) 'default_lambda', (2) 'lambda_from_grid_final_values',
# (3) 'lambda_from_profile_gravitational',
# (4) 'lambda_from_profile_gravitational_plus_internal',
# (5) 'lambda_from_profile_gravitational_plus_internal_minus_recombination'
common_envelope_lambda_default=0.5, # float in (0, inf) used only for option (1)
common_envelope_option_for_HG_star="optimistic", # 'optimistic', 'pessimistic'
common_envelope_alpha_thermal=1.0, # float in (0, inf) used only for option for (4), (5)
core_definition_H_fraction=0.1, # 0.01, 0.1, 0.3
core_definition_He_fraction=0.1, # 0.1
CEE_tolerance_err=0.001, # float (0, inf)
common_envelope_option_after_succ_CEE = 'core_not_replaced_noMT', # 'core_not_replaced_noMT' 'core_replaced_noMT' 'core_not_replaced_stableMT' 'core_not_replaced_windloss'
verbose = False, # True False
)
SN_STEP = dict(
mechanism='Patton&Sukhbold20-engine', # 'direct', Fryer+12-rapid', 'Fryer+12-delayed', 'Sukhbold+16-engine', 'Patton&Sukhbold20-engine'
engine='N20', # 'N20' for 'Sukhbold+16-engine', 'Patton&Sukhbold20-engine' or None for the others
PISN="Marchant+19", # None, "Marchant+19"
ECSN="Podsiadlowksi+04", # "Tauris+15", "Podsiadlowksi+04"
conserve_hydrogen_envelope=True,
max_neutrino_mass_loss=0.5, # float (0,inf)
kick=True, # True, False
kick_normalisation='one_over_mass', # "one_minus_fallback", "one_over_mass", "NS_one_minus_fallback_BH_one", "one", "zero"
sigma_kick_CCSN_NS=265.0, # float (0,inf)
sigma_kick_CCSN_BH=265.0, # float (0,inf)
sigma_kick_ECSN=20.0, # float (0,inf)
max_NS_mass=2.5, # float (0,inf)
use_interp_values=True, # True, False
use_profiles=True, # True, False
use_core_masses=True, # True, False
approx_at_he_depletion=False, # True, False
verbose = False, # True False
)
DCO_STEP = dict(
n_o_steps_interval = None,
)
END_STEP = {}
# FLOW CHART CONFIGURATION
sim_kwargs = dict(
flow = (flow_chart, {}),
step_HMS_HMS = (MS_MS_step, MESA_STEP),
step_CO_HeMS = (CO_HeMS_step, MESA_STEP),
step_CO_HMS_RLO = (CO_HMS_RLO_step, MESA_STEP),
step_CO_HeMS_RLO = (CO_HeMS_RLO_step, MESA_STEP),
step_detached = (detached_step, DETACHED_STEP),
step_disrupted = (DisruptedStep, DISRUPTED_STEP),
step_merged = (MergedStep, {}),
step_CE = (StepCEE, CE_STEP),
step_SN = (StepSN, SN_STEP),
step_dco = (DoubleCO, DCO_STEP),
step_end = (step_end, END_STEP),
extra_hooks = [(StepNamesHooks, {})]
)
sim_prop = SimulationProperties(**sim_kwargs)
# SIMULATION CONFIGURATION
kwargs = dict(
file_path='./batches/',
optimize_ram=False,
ram_per_cpu=3., # limit ram usage at 3GB
dump_rate=1000, # limit batch size
metallicity = 0.0001, # 1e+00, 1e-01, 1e-02, 1e-03, 1e-04
number_of_binaries=10, # int
star_formation='burst', # 'constant' 'burst' 'custom_linear' 'custom_log10' 'custom_linear_histogram' 'custom_log10_histogram'
max_simulation_time=13.8e9, # float (0,inf)
primary_mass_scheme='Kroupa2001', # 'Salpeter', 'Kroupa1993', 'Kroupa2001'
primary_mass_min=6.5, # float (0,130)
primary_mass_max=250., # float (0,130)
secondary_mass_scheme='flat_mass_ratio', # 'flat_mass_ratio', 'q=1'
secondary_mass_min=0.35, # float (0,130)
secondary_mass_max=250., # float (0,130)
orbital_scheme = 'period', # 'separation', 'period'
orbital_period_scheme = 'Sana+12_period_extended', # used only for orbital_scheme = 'period'
orbital_period_min = 1., # float (0,inf)
orbital_period_max = 1000., # float (0,inf)
#orbital_separation_scheme='log_uniform', # used only for orbital_scheme = 'separation', 'log_uniform', 'log_normal'
#orbital_separation_min=5., # float (0,inf)
#orbital_separation_max=1e5, # float (0,inf)
#log_orbital_separation_mean=None, # float (0,inf) used only for orbital_separation_scheme ='log_normal'
#log_orbital_separation_sigma=None, # float (0,inf) used only for orbital_separation_scheme ='log_normal'
eccentricity_sche='zero', # 'zero' 'thermal' 'uniform'
# IMPORT CUSTOM HOOKS
extra_columns={'step_names':'string'}, # 'step_times' with from posydon.binary_evol.simulationproperties import TimingHooks
# LIST BINARY PROPERTIES TO SAVE
only_select_columns=[
'state',
'event',
'time',
#'separation',
'orbital_period',
'eccentricity',
#'V_sys',
#'rl_relative_overflow_1',
#'rl_relative_overflow_2',
'lg_mtransfer_rate',
#'mass_transfer_case',
#'trap_radius',
#'acc_radius',
#'t_sync_rad_1',
#'t_sync_conv_1',
#'t_sync_rad_2',
#'t_sync_conv_2',
#'nearest_neighbour_distance',
],
scalar_names=[
'interp_class_HMS_HMS',
'interp_class_CO_HMS_RLO',
'interp_class_CO_HeMS',
'interp_class_CO_HeMS_RLO'
],
# LIST STAR PROPERTIES TO SAVE
include_S1=True , # True, False
S1_kwargs=dict(only_select_columns=[
'state',
'metallicity',
'mass',
'log_R',
'log_L',
'lg_mdot',
#'lg_system_mdot',
#'lg_wind_mdot',
'he_core_mass',
'he_core_radius',
#'c_core_mass',
#'c_core_radius',
#'o_core_mass',
#'o_core_radius',
'co_core_mass',
'co_core_radius',
'center_h1',
'center_he4',
#'center_c12',
#'center_n14',
#'center_o16',
'surface_h1',
'surface_he4',
#'surface_c12',
#'surface_n14',
#'surface_o16',
#'log_LH',
#'log_LHe',
#'log_LZ',
#'log_Lnuc',
#'c12_c12',
#'center_gamma',
#'avg_c_in_c_core',
#'surf_avg_omega',
'surf_avg_omega_div_omega_crit',
#'total_moment_of_inertia',
#'log_total_angular_momentum',
'spin',
#'conv_env_top_mass',
#'conv_env_bot_mass',
#'conv_env_top_radius',
#'conv_env_bot_radius',
#'conv_env_turnover_time_g',
#'conv_env_turnover_time_l_b',
#'conv_env_turnover_time_l_t',
#'envelope_binding_energy',
#'mass_conv_reg_fortides',
#'thickness_conv_reg_fortides',
#'radius_conv_reg_fortides',
#'lambda_CE_1cent',
#'lambda_CE_10cent',
#'lambda_CE_30cent',
#'lambda_CE_pure_He_star_10cent',
#'profile',
],
scalar_names=['natal_kick_array',
'SN_type',
#'f_fb',
#'spin_orbit_tilt',
]),
# LIST STAR PROPERTIES TO SAVE
include_S2=True, # True, False
S2_kwargs=dict(only_select_columns=[
'state',
'metallicity',
'mass',
'log_R',
'log_L',
'lg_mdot',
#'lg_system_mdot',
#'lg_wind_mdot',
'he_core_mass',
'he_core_radius',
#'c_core_mass',
#'c_core_radius',
#'o_core_mass',
#'o_core_radius',
'co_core_mass',
'co_core_radius',
'center_h1',
'center_he4',
#'center_c12',
#'center_n14',
#'center_o16',
'surface_h1',
'surface_he4',
#'surface_c12',
#'surface_n14',
#'surface_o16',
#'log_LH',
#'log_LHe',
#'log_LZ',
#'log_Lnuc',
#'c12_c12',
#'center_gamma',
#'avg_c_in_c_core',
#'surf_avg_omega',
'surf_avg_omega_div_omega_crit',
#'total_moment_of_inertia',
#'log_total_angular_momentum',
'spin',
#'conv_env_top_mass',
#'conv_env_bot_mass',
#'conv_env_top_radius',
#'conv_env_bot_radius',
#'conv_env_turnover_time_g',
#'conv_env_turnover_time_l_b',
#'conv_env_turnover_time_l_t',
#'envelope_binding_energy',
#'mass_conv_reg_fortides',
#'thickness_conv_reg_fortides',
#'radius_conv_reg_fortides',
#'lambda_CE_1cent',
#'lambda_CE_10cent',
#'lambda_CE_30cent',
#'lambda_CE_pure_He_star_10cent',
#'profile',
],
scalar_names=['natal_kick_array',
'SN_type',
#'f_fb',
#'spin_orbit_tilt',
]),
)
def run_simulation(sim_prop, kwargs, file=None, indices=None):
# create binaries
pop = BinaryPopulation(entropy=None,
population_properties=sim_prop,
file_name=file,
**kwargs)
sim_prop.load_steps(verbose=True)
# evolve binaries
if file is not None:
kwargs['from_hdf'] = True
kwargs['indices'] = indices
pop.evolve(breakdown_to_df=False, tqdm=True, **kwargs)
# save binaries
try:
pop.save('./population.h5', **kwargs)
except:
pass
return pop
if __name__ == '__main__' :
pop = run_simulation(sim_prop, kwargs)
Letβs say we want to debug binary_index that follows.
[ ]:
i = 0
Re-evolving a Binary Starο
[ ]:
col = ['time', 'step_names', 'state', 'event', 'orbital_period', 'eccentricity', 'S1_state', 'S2_state', 'S1_mass', 'S2_mass']
pop[i].to_df(extra_columns={'step_names':'string'})[col]
[ ]:
pop[i].restore()
[ ]:
pop[i].to_df()
Even thought we reset the binary, the natal kick properties stored in the natal kick array are not reset. This is because we want to be able to re-evolve the binary to the same final state. In case you want to reset the natal kick array values, then set the list to None values.
[ ]:
# kick magnitude km/s, azimuthal angle rad, polar angle rad and mean anomaly rad (TODO check documentation)
print(pop[i].star_1.natal_kick_array)
print(pop[i].star_2.natal_kick_array)
[ ]:
pop[i].evolve()
[ ]:
pop[i].to_df(extra_columns={'step_names':'string'})[col]
Re-evolving a Binary Star with a different kickο
[ ]:
pop[i].restore()
pop[i].to_df()
[ ]:
pop[i].star_1.natal_kick_array = [1000., 0.858129274334538, 1.9157148786534735, 1.8675467897282945]
pop[i].star_2.natal_kick_array = [None, None, None, None] # default
print(pop[i].star_1.natal_kick_array)
print(pop[i].star_2.natal_kick_array)
[ ]:
pop[i].evolve()
pop[i].to_df(extra_columns={'step_names':'string'})[col]
Debugging a Binary Star from a population h5 fileο
Assume we run the population somewhere else and you have the h5 file associated with it and you want to load the population in memory and re-evolve a binary. Here is how you do it.
[ ]:
del pop # delete the above population
[ ]:
# let's load up and rerun the binary index i
pop = run_simulation(sim_prop, kwargs, file='./population.h5', indices=[i])
[ ]:
pop[0].to_df(extra_columns={'step_names':'string'})[col]
[ ]:
# do the same thing as above
pop[0].restore()
print(pop[i].star_1.natal_kick_array)
print(pop[i].star_2.natal_kick_array)
pop[i].star_1.natal_kick_array = [1000., 0.858129274334538, 1.9157148786534735, 1.8675467897282945]
pop[i].star_2.natal_kick_array = [None, None, None, None] # default
print(pop[i].star_1.natal_kick_array)
print(pop[i].star_2.natal_kick_array)
pop[i].evolve()
pop[i].to_df(extra_columns={'step_names':'string'})[col]
Evolving a Binary Star starting from an arbitrary stateο
If you want to evolve a custom binary star, you can do so by crafting yourself the BinaryStar object and providing the SimulationProperties object. For example, letβs evolve a neutron star with a low mass helium star in roche lobe overflow.
Caution!
While you can evolve a binary from an arbitrary state, you will need to provide data on the internal structure of the star, if youβre starting in the detached step. Otherwise the matching to a single star model will not work!
[ ]:
from posydon.binary_evol.singlestar import SingleStar
from posydon.binary_evol.binarystar import BinaryStar
sim_prop.load_steps(verbose=True) # load steps
[ ]:
binary = BinaryStar(star_1 = SingleStar(**{'state' : 'NS',
'mass' : 1.1,
'spin' : 0.,}),
star_2 = SingleStar(**{'state' : 'stripped_He_Core_He_burning',
'mass' : 2.5,
'natal_kick_array' : [10., 0., 0., 0.]}),
**{'time' : 0.,
'state' : 'RLO2',
'event' : 'oRLO2',
'orbital_period' : 1.,
'eccentricity' : 0.},
properties = sim_prop,
)
binary.evolve()
binary.to_df(extra_columns={'step_names':'string'})[col]
Mapping a MESA step to the detached_stepο
You might want to use the detached_step instetad of a MESA step. Here is how you do it. As an example we replace the HMS-HMS MESA step with the detached step. Notice that the detached_step does not evolve the binary system post onset of RLO1. Hence, unless the flow_chart support such binary event, binary state and stellar states, else these systems will not be furthter evovlved.
[ ]:
# MODIFIED FLOW CHART CONFIGURATION
sim_kwargs = dict(
flow = (flow_chart, {}),
step_HMS_HMS = (detached_step, DETACHED_STEP), # <-- map this step to detached
step_CO_HeMS = (CO_HeMS_step, MESA_STEP),
step_CO_HMS_RLO = (CO_HMS_RLO_step, MESA_STEP),
step_CO_HeMS_RLO = (CO_HeMS_RLO_step, MESA_STEP),
step_detached = (detached_step, DETACHED_STEP),
step_disrupted = (DisruptedStep, DISRUPTED_STEP),
step_merged = (MergedStep, {}),
step_CE = (StepCEE, CE_STEP),
step_SN = (StepSN, SN_STEP),
step_dco = (DoubleCO, DCO_STEP),
step_end = (step_end, END_STEP),
extra_hooks = [(StepNamesHooks, {})]
)
sim_prop = SimulationProperties(**sim_kwargs)
kwargs['from_hdf'] = False
if __name__ == '__main__' :
pop = run_simulation(sim_prop, kwargs)
[ ]:
i = 0
pop[i].to_df(extra_columns={'step_names':'string'})[col]
Crating your own stepο
Here, you learn how to create your own step. As an example we create a step that does a fake common envelope by taking half the orbital period.
[ ]:
class my_CE_step(object):
"""Compute a fake CE event."""
def __init__(self, verbose=False):
self.verbose = verbose
def __call__(self, binary):
if self.verbose:
print('The orbital separation post CE is half the pre CE orbital separation!')
# Determine which star is the donor and which is the companion
if binary.event in ["oCE1", "oDoubleCE1"]:
donor_star = binary.star_1
comp_star = binary.star_2
elif binary.event in ["oCE2", "oDoubleCE2"]:
donor_star = binary.star_2
comp_star = binary.star_1
else:
raise ValueError("CEE does not apply if `event` is not "
"`oCE1`, 'oDoubleCE1' or `oCE2`, 'oDoubleCE1'")
binary.orbital_period /= 2.
donor_star.mass = donor_star.he_core_mass # lose envelope
donor_star.state = donor_star.state.replace('H-rich', 'stripped_He')
binary.state = 'detached'
binary.event = None
# MODIFIED FLOW CHART CONFIGURATION
sim_kwargs = dict(
flow = (flow_chart, {}),
step_HMS_HMS = (MS_MS_step, MESA_STEP),
step_CO_HeMS = (CO_HeMS_step, MESA_STEP),
step_CO_HMS_RLO = (CO_HMS_RLO_step, MESA_STEP),
step_CO_HeMS_RLO = (CO_HeMS_RLO_step, MESA_STEP),
step_detached = (detached_step, DETACHED_STEP),
step_disrupted = (DisruptedStep, DISRUPTED_STEP),
step_merged = (MergedStep, {}),
step_CE = (my_CE_step, {}), # <-- map this step to my fake CE step
step_SN = (StepSN, SN_STEP),
step_dco = (DoubleCO, DCO_STEP),
step_end = (step_end, END_STEP),
extra_hooks = [(StepNamesHooks, {})]
)
sim_prop = SimulationProperties(**sim_kwargs)
pop = run_simulation(sim_prop, kwargs)
[ ]:
i = 0
pop[i].to_df(extra_columns={'step_names':'string'})[col]