{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Custom POSYDON steps and flow chart\n", "\n", "The evolution of binaries in POSYDON is determined by the flow chart.\n", "It uses the stellar states, `S1_state` and `S2_state`, the binary `state` and `event` to determine to which step the binary goes next.\n", "\n", "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`.\n", "\n", "If your specific need requires, it is possible to change the flow and introduce your own steps.\n", "This is possible through\n", "\n", "1. Importing an additional flow chart and step from the `population_params.ini` file\n", "2. Changing a step inside a notebook\n", "\n", "\n", "## Population_params.ini changes\n", "\n", "If you want to run large populations with an adapted flowchart or custom step, this is the location to change the loaded steps.\n", "\n", "### Importing an adapted flowchart\n", "\n", "One of the first things the `population_params.ini` file defines is the flow chart.\n", "By default, the POSYDON default `flow_chart` is imported:\n", "```\n", "[flow]\n", " import = ['posydon.binary_evol.flow_chart', 'flow_chart']\n", " # builtin posydon flow\n", " absolute_import = None\n", " # If given, use an absolute filepath to user defined flow: ['', '']\n", "```\n", "\n", "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.\n", "\n", "Additionally, you can put the python file into the `user_modules` folder to make it part of the POSYDON namespace.\n", "\n", "```\n", "[flow]\n", " import = None\n", " # builtin posydon flow\n", " absolute_import = ['custom_step_and_flow.py', 'end_flow_chart']\n", " # If given, use an absolute filepath to user defined flow: ['', '']\n", "```\n", "\n", "We can check this change by loading the Simulation Property arguments.\n", "Below you can see the example, where we send all the states straight to `step_end`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from posydon.popsyn.io import simprop_kwargs_from_ini\n", "from posydon.binary_evol.simulationproperties import SimulationProperties\n", "\n", "sim_kwargs = simprop_kwargs_from_ini('population_params.ini')\n", "sim_prop = SimulationProperties(**sim_kwargs)\n", "sim_prop.flow[0]()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim_kwargs = simprop_kwargs_from_ini('population_params_end.ini')\n", "sim_prop = SimulationProperties(**sim_kwargs)\n", "sim_prop.flow[0]()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Changing an existing step to a custom step\n", "\n", "Similarly, you can change one of standard POSYDON steps to a custom one.\n", "We will replace the common envelope step by one that halves the orbital period.\n", "\n", "```\n", "[step_CE]\n", " import = ['posydon.binary_evol.CE.step_CEE', 'StepCEE']\n", " # builtin posydon step\n", " absolute_import = None\n", " # If given, use an absolute filepath to user defined step: ['', '']\n", "```\n", "\n", "```\n", "[step_CE]\n", " import = None \n", " # builtin posydon step\n", " absolute_import = ['custom_step_and_flow.py', 'my_CE_step']\n", " # If given, use an absolute filepath to user defined step: ['', '']\n", "```\n", "\n", "You can see the step class in the example file [custom_step_and_flow.py](custom_step_and_flow.py)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim_kwargs = simprop_kwargs_from_ini('population_params_end.ini')\n", "sim_kwargs['step_CE']\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adapting the flow and steps inside your notebook\n", "\n", "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.\n", "\n", "This is definitely not the easiest way to change option within POSYDON." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from posydon.popsyn.binarypopulation import BinaryPopulation\n", "from posydon.binary_evol.simulationproperties import SimulationProperties\n", "from posydon.binary_evol.flow_chart import flow_chart\n", "from posydon.binary_evol.MESA.step_mesa import CO_HeMS_step, MS_MS_step, CO_HMS_RLO_step, CO_HeMS_RLO_step\n", "from posydon.binary_evol.DT.step_detached import detached_step\n", "from posydon.binary_evol.DT.step_disrupted import DisruptedStep\n", "from posydon.binary_evol.DT.step_merged import MergedStep\n", "from posydon.binary_evol.CE.step_CEE import StepCEE\n", "from posydon.binary_evol.SN.step_SN import StepSN\n", "from posydon.binary_evol.DT.double_CO import DoubleCO\n", "from posydon.binary_evol.step_end import step_end\n", "from posydon.binary_evol.simulationproperties import StepNamesHooks\n", "\n", "# STEP CUSTOMISATION\n", "MESA_STEP = dict(\n", " interpolation_path = None, # found by default\n", " interpolation_filename = None, # found by default\n", " interpolation_method = 'nearest_neighbour', # 'nearest_neighbour' 'linear3c_kNN' '1NN_1NN'\n", " save_initial_conditions = True, # only for interpolation_method='nearest_neighbour'\n", " track_interpolation = False, # True False\n", " stop_method = 'stop_at_max_time', # 'stop_at_end' 'stop_at_max_time' 'stop_at_condition'\n", " stop_star = 'star_1', # only for stop_method='stop_at_condition' 'star_1' 'star_2'\n", " stop_var_name = None, # only for stop_method='stop_at_condition' str\n", " stop_value = None, # only for stop_method='stop_at_condition' float\n", " stop_interpolate = True, # True False\n", " verbose = False, # True False\n", ")\n", "\n", "DETACHED_STEP = dict(\n", " matching_method = 'minimize', #'minimize' 'root'\n", " do_wind_loss = True, # True False\n", " do_tides = True, # True False\n", " do_gravitational_radiation = True, # True False\n", " do_magnetic_braking = True, # True False\n", " do_stellar_evolution_and_spin_from_winds = True, # True False\n", " RLO_orbit_at_orbit_with_same_am = False, # True False\n", " verbose = False, # True False\n", ")\n", "\n", "DISRUPTED_STEP = dict(\n", " grid_name_Hrich=None,\n", " grid_name_strippedHe=None,\n", " metallicity=None,\n", " dt=None,\n", " n_o_steps_history=None,\n", " matching_method=\"minimize\",\n", " initial_mass=None,\n", " rootm=None,\n", " verbose=False,\n", " do_wind_loss=True,\n", " do_tides=True,\n", " do_gravitational_radiation=True,\n", " do_magnetic_braking=True,\n", " magnetic_braking_mode=\"RVJ83\",\n", " do_stellar_evolution_and_spin_from_winds=True,\n", " RLO_orbit_at_orbit_with_same_am=False,\n", " list_for_matching_HMS=None,\n", " list_for_matching_postMS=None,\n", " list_for_matching_HeStar=None\n", ")\n", "\n", "CE_STEP = dict(\n", " prescription='alpha-lambda', # 'alpha-lambda'\n", " common_envelope_efficiency=1.0, # float in (0, inf)\n", " common_envelope_option_for_lambda='lambda_from_grid_final_values', # (1) 'default_lambda', (2) 'lambda_from_grid_final_values',\n", " # (3) 'lambda_from_profile_gravitational',\n", " # (4) 'lambda_from_profile_gravitational_plus_internal',\n", " # (5) 'lambda_from_profile_gravitational_plus_internal_minus_recombination'\n", " common_envelope_lambda_default=0.5, # float in (0, inf) used only for option (1)\n", " common_envelope_option_for_HG_star=\"optimistic\", # 'optimistic', 'pessimistic'\n", " common_envelope_alpha_thermal=1.0, # float in (0, inf) used only for option for (4), (5)\n", " core_definition_H_fraction=0.1, # 0.01, 0.1, 0.3\n", " core_definition_He_fraction=0.1, # 0.1\n", " CEE_tolerance_err=0.001, # float (0, inf)\n", " 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'\n", " verbose = False, # True False\n", ")\n", "\n", "SN_STEP = dict(\n", " mechanism='Patton&Sukhbold20-engine', # 'direct', Fryer+12-rapid', 'Fryer+12-delayed', 'Sukhbold+16-engine', 'Patton&Sukhbold20-engine'\n", " engine='N20', # 'N20' for 'Sukhbold+16-engine', 'Patton&Sukhbold20-engine' or None for the others\n", " PISN=\"Marchant+19\", # None, \"Marchant+19\"\n", " ECSN=\"Podsiadlowksi+04\", # \"Tauris+15\", \"Podsiadlowksi+04\"\n", " conserve_hydrogen_envelope=True,\n", " max_neutrino_mass_loss=0.5, # float (0,inf)\n", " kick=True, # True, False\n", " kick_normalisation='one_over_mass', # \"one_minus_fallback\", \"one_over_mass\", \"NS_one_minus_fallback_BH_one\", \"one\", \"zero\"\n", " sigma_kick_CCSN_NS=265.0, # float (0,inf)\n", " sigma_kick_CCSN_BH=265.0, # float (0,inf)\n", " sigma_kick_ECSN=20.0, # float (0,inf)\n", " max_NS_mass=2.5, # float (0,inf)\n", " use_interp_values=True, # True, False\n", " use_profiles=True, # True, False\n", " use_core_masses=True, # True, False\n", " approx_at_he_depletion=False, # True, False\n", " verbose = False, # True False\n", ")\n", "\n", "DCO_STEP = dict(\n", " n_o_steps_interval = None,\n", ")\n", "\n", "END_STEP = {}\n", "\n", "# FLOW CHART CONFIGURATION\n", "sim_kwargs = dict(\n", " flow = (flow_chart, {}),\n", " step_HMS_HMS = (MS_MS_step, MESA_STEP),\n", " step_CO_HeMS = (CO_HeMS_step, MESA_STEP),\n", " step_CO_HMS_RLO = (CO_HMS_RLO_step, MESA_STEP),\n", " step_CO_HeMS_RLO = (CO_HeMS_RLO_step, MESA_STEP),\n", " step_detached = (detached_step, DETACHED_STEP),\n", " step_disrupted = (DisruptedStep, DISRUPTED_STEP),\n", " step_merged = (MergedStep, {}),\n", " step_CE = (StepCEE, CE_STEP),\n", " step_SN = (StepSN, SN_STEP),\n", " step_dco = (DoubleCO, DCO_STEP),\n", " step_end = (step_end, END_STEP),\n", " extra_hooks = [(StepNamesHooks, {})]\n", ")\n", "\n", "sim_prop = SimulationProperties(**sim_kwargs)\n", "\n", "# SIMULATION CONFIGURATION\n", "kwargs = dict(\n", " file_path='./batches/',\n", " optimize_ram=False,\n", " ram_per_cpu=3., # limit ram usage at 3GB\n", " dump_rate=1000, # limit batch size\n", "\n", " metallicity = 0.0001, # 1e+00, 1e-01, 1e-02, 1e-03, 1e-04\n", " number_of_binaries=10, # int\n", " star_formation='burst', # 'constant' 'burst' 'custom_linear' 'custom_log10' 'custom_linear_histogram' 'custom_log10_histogram'\n", " max_simulation_time=13.8e9, # float (0,inf)\n", "\n", " primary_mass_scheme='Kroupa2001', # 'Salpeter', 'Kroupa1993', 'Kroupa2001'\n", " primary_mass_min=6.5, # float (0,130)\n", " primary_mass_max=250., # float (0,130)\n", " secondary_mass_scheme='flat_mass_ratio', # 'flat_mass_ratio', 'q=1'\n", " secondary_mass_min=0.35, # float (0,130)\n", " secondary_mass_max=250., # float (0,130)\n", " orbital_scheme = 'period', # 'separation', 'period'\n", " orbital_period_scheme = 'Sana+12_period_extended', # used only for orbital_scheme = 'period'\n", " orbital_period_min = 1., # float (0,inf)\n", " orbital_period_max = 1000., # float (0,inf)\n", " #orbital_separation_scheme='log_uniform', # used only for orbital_scheme = 'separation', 'log_uniform', 'log_normal'\n", " #orbital_separation_min=5., # float (0,inf)\n", " #orbital_separation_max=1e5, # float (0,inf)\n", " #log_orbital_separation_mean=None, # float (0,inf) used only for orbital_separation_scheme ='log_normal'\n", " #log_orbital_separation_sigma=None, # float (0,inf) used only for orbital_separation_scheme ='log_normal'\n", " eccentricity_sche='zero', # 'zero' 'thermal' 'uniform'\n", "\n", " # IMPORT CUSTOM HOOKS\n", " extra_columns={'step_names':'string'}, # 'step_times' with from posydon.binary_evol.simulationproperties import TimingHooks\n", "\n", " # LIST BINARY PROPERTIES TO SAVE\n", " only_select_columns=[\n", " 'state',\n", " 'event',\n", " 'time',\n", " #'separation',\n", " 'orbital_period',\n", " 'eccentricity',\n", " #'V_sys',\n", " #'rl_relative_overflow_1',\n", " #'rl_relative_overflow_2',\n", " 'lg_mtransfer_rate',\n", " #'mass_transfer_case',\n", " #'trap_radius',\n", " #'acc_radius',\n", " #'t_sync_rad_1',\n", " #'t_sync_conv_1',\n", " #'t_sync_rad_2',\n", " #'t_sync_conv_2',\n", " #'nearest_neighbour_distance',\n", " ],\n", " scalar_names=[\n", " 'interp_class_HMS_HMS',\n", " 'interp_class_CO_HMS_RLO',\n", " 'interp_class_CO_HeMS',\n", " 'interp_class_CO_HeMS_RLO'\n", " ],\n", "\n", " # LIST STAR PROPERTIES TO SAVE\n", " include_S1=True , # True, False\n", " S1_kwargs=dict(only_select_columns=[\n", " 'state',\n", " 'metallicity',\n", " 'mass',\n", " 'log_R',\n", " 'log_L',\n", " 'lg_mdot',\n", " #'lg_system_mdot',\n", " #'lg_wind_mdot',\n", " 'he_core_mass',\n", " 'he_core_radius',\n", " #'c_core_mass',\n", " #'c_core_radius',\n", " #'o_core_mass',\n", " #'o_core_radius',\n", " 'co_core_mass',\n", " 'co_core_radius',\n", " 'center_h1',\n", " 'center_he4',\n", " #'center_c12',\n", " #'center_n14',\n", " #'center_o16',\n", " 'surface_h1',\n", " 'surface_he4',\n", " #'surface_c12',\n", " #'surface_n14',\n", " #'surface_o16',\n", " #'log_LH',\n", " #'log_LHe',\n", " #'log_LZ',\n", " #'log_Lnuc',\n", " #'c12_c12',\n", " #'center_gamma',\n", " #'avg_c_in_c_core',\n", " #'surf_avg_omega',\n", " 'surf_avg_omega_div_omega_crit',\n", " #'total_moment_of_inertia',\n", " #'log_total_angular_momentum',\n", " 'spin',\n", " #'conv_env_top_mass',\n", " #'conv_env_bot_mass',\n", " #'conv_env_top_radius',\n", " #'conv_env_bot_radius',\n", " #'conv_env_turnover_time_g',\n", " #'conv_env_turnover_time_l_b',\n", " #'conv_env_turnover_time_l_t',\n", " #'envelope_binding_energy',\n", " #'mass_conv_reg_fortides',\n", " #'thickness_conv_reg_fortides',\n", " #'radius_conv_reg_fortides',\n", " #'lambda_CE_1cent',\n", " #'lambda_CE_10cent',\n", " #'lambda_CE_30cent',\n", " #'lambda_CE_pure_He_star_10cent',\n", " #'profile',\n", " ],\n", " scalar_names=['natal_kick_array',\n", " 'SN_type',\n", " #'f_fb',\n", " #'spin_orbit_tilt',\n", " ]),\n", "\n", " # LIST STAR PROPERTIES TO SAVE\n", " include_S2=True, # True, False\n", " S2_kwargs=dict(only_select_columns=[\n", " 'state',\n", " 'metallicity',\n", " 'mass',\n", " 'log_R',\n", " 'log_L',\n", " 'lg_mdot',\n", " #'lg_system_mdot',\n", " #'lg_wind_mdot',\n", " 'he_core_mass',\n", " 'he_core_radius',\n", " #'c_core_mass',\n", " #'c_core_radius',\n", " #'o_core_mass',\n", " #'o_core_radius',\n", " 'co_core_mass',\n", " 'co_core_radius',\n", " 'center_h1',\n", " 'center_he4',\n", " #'center_c12',\n", " #'center_n14',\n", " #'center_o16',\n", " 'surface_h1',\n", " 'surface_he4',\n", " #'surface_c12',\n", " #'surface_n14',\n", " #'surface_o16',\n", " #'log_LH',\n", " #'log_LHe',\n", " #'log_LZ',\n", " #'log_Lnuc',\n", " #'c12_c12',\n", " #'center_gamma',\n", " #'avg_c_in_c_core',\n", " #'surf_avg_omega',\n", " 'surf_avg_omega_div_omega_crit',\n", " #'total_moment_of_inertia',\n", " #'log_total_angular_momentum',\n", " 'spin',\n", " #'conv_env_top_mass',\n", " #'conv_env_bot_mass',\n", " #'conv_env_top_radius',\n", " #'conv_env_bot_radius',\n", " #'conv_env_turnover_time_g',\n", " #'conv_env_turnover_time_l_b',\n", " #'conv_env_turnover_time_l_t',\n", " #'envelope_binding_energy',\n", " #'mass_conv_reg_fortides',\n", " #'thickness_conv_reg_fortides',\n", " #'radius_conv_reg_fortides',\n", " #'lambda_CE_1cent',\n", " #'lambda_CE_10cent',\n", " #'lambda_CE_30cent',\n", " #'lambda_CE_pure_He_star_10cent',\n", " #'profile',\n", " ],\n", " scalar_names=['natal_kick_array',\n", " 'SN_type',\n", " #'f_fb',\n", " #'spin_orbit_tilt',\n", " ]),\n", ")\n", "\n", "\n", "def run_simulation(sim_prop, kwargs, file=None, indices=None):\n", "\n", "\n", " # create binaries\n", " pop = BinaryPopulation(entropy=None,\n", " population_properties=sim_prop,\n", " file_name=file,\n", " **kwargs)\n", " \n", " sim_prop.load_steps(verbose=True)\n", "\n", " # evolve binaries\n", " if file is not None:\n", " kwargs['from_hdf'] = True\n", " kwargs['indices'] = indices\n", "\n", " pop.evolve(breakdown_to_df=False, tqdm=True, **kwargs)\n", "\n", " # save binaries\n", " try:\n", " pop.save('./population.h5', **kwargs)\n", " except:\n", " pass\n", "\n", " return pop\n", "\n", " \n", "col = ['time', 'step_names', 'state', 'event', 'orbital_period', 'eccentricity', 'S1_state', 'S2_state', 'S1_mass', 'S2_mass']\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adapted flow chart\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# MODIFIED FLOW CHART CONFIGURATION\n", "sim_kwargs = dict(\n", " flow = (flow_chart, {}),\n", " step_HMS_HMS = (detached_step, DETACHED_STEP), # <-- map this step to detached\n", " step_CO_HeMS = (CO_HeMS_step, MESA_STEP),\n", " step_CO_HMS_RLO = (CO_HMS_RLO_step, MESA_STEP),\n", " step_CO_HeMS_RLO = (CO_HeMS_RLO_step, MESA_STEP),\n", " step_detached = (detached_step, DETACHED_STEP),\n", " step_disrupted = (DisruptedStep, DISRUPTED_STEP),\n", " step_merged = (MergedStep, {}),\n", " step_CE = (StepCEE, CE_STEP),\n", " step_SN = (StepSN, SN_STEP),\n", " step_dco = (DoubleCO, DCO_STEP),\n", " step_end = (step_end, END_STEP),\n", " extra_hooks = [(StepNamesHooks, {})]\n", ")\n", "\n", "sim_prop = SimulationProperties(**sim_kwargs)\n", "\n", "kwargs['from_hdf'] = False\n", "\n", "pop = run_simulation(sim_prop, kwargs)\n", "# output the first binary\n", "i = 0\n", "pop[i].to_df(extra_columns={'step_names':'string'})[col]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adapted step\n", "\n", "We implement a custom common envelope step, which is the same as the on in the `custom_step_and_flow.py` custom step." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class my_CE_step(object):\n", " \"\"\"Compute a fake CE event.\"\"\"\n", "\n", " def __init__(self, verbose=False):\n", " self.verbose = verbose\n", "\n", " def __call__(self, binary):\n", "\n", " if self.verbose:\n", " print('The orbital separation post CE is half the pre CE orbital separation!')\n", "\n", " # Determine which star is the donor and which is the companion\n", " if binary.event in [\"oCE1\", \"oDoubleCE1\"]:\n", " donor_star = binary.star_1\n", " comp_star = binary.star_2\n", " elif binary.event in [\"oCE2\", \"oDoubleCE2\"]:\n", " donor_star = binary.star_2\n", " comp_star = binary.star_1\n", " else:\n", " raise ValueError(\"CEE does not apply if `event` is not \"\n", " \"`oCE1`, 'oDoubleCE1' or `oCE2`, 'oDoubleCE1'\")\n", "\n", " binary.orbital_period /= 2.\n", " donor_star.mass = donor_star.he_core_mass # lose envelope\n", " donor_star.state = donor_star.state.replace('H-rich', 'stripped_He')\n", " binary.state = 'detached'\n", " binary.event = None\n", " \n", " \n", "# MODIFIED FLOW CHART CONFIGURATION\n", "sim_kwargs = dict(\n", " flow = (flow_chart, {}),\n", " step_HMS_HMS = (MS_MS_step, MESA_STEP),\n", " step_CO_HeMS = (CO_HeMS_step, MESA_STEP),\n", " step_CO_HMS_RLO = (CO_HMS_RLO_step, MESA_STEP),\n", " step_CO_HeMS_RLO = (CO_HeMS_RLO_step, MESA_STEP),\n", " step_detached = (detached_step, DETACHED_STEP),\n", " step_disrupted = (DisruptedStep, DISRUPTED_STEP),\n", " step_merged = (MergedStep, {}),\n", " step_CE = (my_CE_step, {}), # <-- map this step to my fake CE step\n", " step_SN = (StepSN, SN_STEP),\n", " step_dco = (DoubleCO, DCO_STEP),\n", " step_end = (step_end, END_STEP),\n", " extra_hooks = [(StepNamesHooks, {})]\n", ")\n", "\n", "sim_prop = SimulationProperties(**sim_kwargs)\n", "\n", "pop = run_simulation(sim_prop, kwargs)\n", "i = 0\n", "pop[i].to_df(extra_columns={'step_names':'string'})[col]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "kernelspec": { "display_name": "posydon_env", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 2 }