Custom POSYDON steps and flow chart
The evolution of binaries in POSYDON is determined by the flow chart. It uses the stellar states, S1_state and S2_state, the binary state and event to determine to which step the binary goes next.
For example, in the default flow chart, the state ('H-rich_Core_C_depleted','H-rich_Core_H_burning', 'detached', 'CC1') will be sent to step_SN.
If your specific need requires, it is possible to change the flow and introduce your own steps. This is possible through
Importing an additional flow chart and step from the
population_params.inifileChanging a step inside a notebook
Population_params.ini changes
If you want to run large populations with an adapted flowchart or custom step, this is the location to change the loaded steps.
Importing an adapted flowchart
One of the first things the population_params.ini file defines is the flow chart. By default, the POSYDON default flow_chart is imported:
[flow]
import = ['posydon.binary_evol.flow_chart', 'flow_chart']
# builtin posydon flow
absolute_import = None
# If given, use an absolute filepath to user defined flow: ['<PATH_TO_PY_FILE>', '<NAME>']
We can change this to import our custom flow. In this example, the binary is immediately send to step_end without any evolution. Make sure to change the path to your absolute path of the example file, if you want to run a population with it.
Additionally, you can put the python file into the user_modules folder to make it part of the POSYDON namespace.
[flow]
import = None
# builtin posydon flow
absolute_import = ['custom_step_and_flow.py', 'end_flow_chart']
# If given, use an absolute filepath to user defined flow: ['<PATH_TO_PY_FILE>', '<NAME>']
We can check this change by loading the Simulation Property arguments. Below you can see the example, where we send all the states straight to step_end.
[ ]:
from posydon.popsyn.io import simprop_kwargs_from_ini
from posydon.binary_evol.simulationproperties import SimulationProperties
sim_kwargs = simprop_kwargs_from_ini('population_params.ini')
sim_prop = SimulationProperties(**sim_kwargs)
sim_prop.flow[0]()
[ ]:
sim_kwargs = simprop_kwargs_from_ini('population_params_end.ini')
sim_prop = SimulationProperties(**sim_kwargs)
sim_prop.flow[0]()
Changing an existing step to a custom step
Similarly, you can change one of standard POSYDON steps to a custom one. We will replace the common envelope step by one that halves the orbital period.
[step_CE]
import = ['posydon.binary_evol.CE.step_CEE', 'StepCEE']
# builtin posydon step
absolute_import = None
# If given, use an absolute filepath to user defined step: ['<PATH_TO_PY_FILE>', '<NAME>']
[step_CE]
import = None
# builtin posydon step
absolute_import = ['custom_step_and_flow.py', 'my_CE_step']
# If given, use an absolute filepath to user defined step: ['<PATH_TO_PY_FILE>', '<NAME>']
You can see the step class in the example file custom_step_and_flow.py.
[ ]:
sim_kwargs = simprop_kwargs_from_ini('population_params_end.ini')
sim_kwargs['step_CE']
Adapting the flow and steps inside your notebook
Although you can always load in the population parameters file into your notebook, maybe you want to create a unique flow inside a notebook. For this you need to set up all the steps manually inside the notebook.
This is definitely not the easiest way to change option within POSYDON.
[ ]:
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
col = ['time', 'step_names', 'state', 'event', 'orbital_period', 'eccentricity', 'S1_state', 'S2_state', 'S1_mass', 'S2_mass']
Adapted flow chart
You might want to use the detached_step instead 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 further evolved.
[ ]:
# 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
pop = run_simulation(sim_prop, kwargs)
# output the first binary
i = 0
pop[i].to_df(extra_columns={'step_names':'string'})[col]
Adapted step
We implement a custom common envelope step, which is the same as the on in the custom_step_and_flow.py custom step.
[ ]:
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]