From 386662c192b38a227cedf5106015e3ddb27034ac Mon Sep 17 00:00:00 2001 From: ezapartas Date: Fri, 24 May 2024 14:22:36 +0300 Subject: [PATCH 1/4] improved mostly step merged with merger object (especially WD mergers), merger center value calcualtion and allowing co_core_mass in detached --- posydon/binary_evol/DT/step_detached.py | 17 ++-- posydon/binary_evol/DT/step_merged.py | 105 ++++++++++++++++-------- posydon/binary_evol/MESA/step_mesa.py | 34 ++++---- 3 files changed, 98 insertions(+), 58 deletions(-) diff --git a/posydon/binary_evol/DT/step_detached.py b/posydon/binary_evol/DT/step_detached.py index a87f4df7d4..e2aa304146 100644 --- a/posydon/binary_evol/DT/step_detached.py +++ b/posydon/binary_evol/DT/step_detached.py @@ -387,7 +387,8 @@ def __init__( "surface_he4", "surface_h1", "log_R", - "center_c12" + "center_c12", + "co_core_mass" ] ) @@ -899,14 +900,15 @@ def sq_diff_function(x): f'{initials[0]:.3f} [Msun],', f'{initials[1]/1e6:.3f} [Myrs]', "\n", "with m(t0), log10(R(t0), center_he(t0), surface_he4(t0), " - "surface_h1(t0), he_core_mass(t0), center_c12(t0) = \n", + "surface_h1(t0), he_core_mass(t0), center_c12(t0), co_core_mass(t0) = \n", f'{self.get_track_val("mass", htrack, *sol.x):.3f}', f'{self.get_track_val("log_R", htrack, *sol.x):.3f}', f'{self.get_track_val("center_he4", htrack, *sol.x):.4f}', f'{self.get_track_val("surface_he4", htrack, *sol.x):.4f}', f'{self.get_track_val("surface_h1", htrack, *sol.x):.4f}', f'{self.get_track_val("he_core_mass", htrack, *sol.x):.3f}', - f'{self.get_track_val("center_c12", htrack, *sol.x):.4f}\n', + f'{self.get_track_val("center_c12", htrack, *sol.x):.4f}', + f'{self.get_track_val("co_core_mass", htrack, *sol.x):.3f}\n', "The same values of the secondary at the end of the previous " "step was = \n", f'{star.mass:.3f}', @@ -915,7 +917,8 @@ def sq_diff_function(x): f'{star.surface_he4:.4f}', f'{star.surface_h1:.4f}', f'{star.he_core_mass:.3f}', - f'{star.center_c12:.4f}' + f'{star.center_c12:.4f}', + f'{star.co_core_mass:.4f}' ) return initials[0], initials[1], htrack @@ -1082,15 +1085,15 @@ def get_star_data(binary, star1, star2, htrack, if self.verbose or self.verbose == 1: print("Matching duration: " f"{t_after_matching-t_before_matching:.6g}") - + if pd.isna(m0) or pd.isna(t0): return None, None, None - + if htrack: self.grid = self.grid_Hrich else: self.grid = self.grid_strippedHe - + # check if m0 is in the grid if m0 < self.grid.grid_mass.min() or m0 > self.grid.grid_mass.max(): set_binary_to_failed(binary) diff --git a/posydon/binary_evol/DT/step_merged.py b/posydon/binary_evol/DT/step_merged.py index 14c5b8478e..ea4336023f 100644 --- a/posydon/binary_evol/DT/step_merged.py +++ b/posydon/binary_evol/DT/step_merged.py @@ -15,6 +15,7 @@ from posydon.utils.common_functions import check_state_of_star from posydon.binary_evol.DT.step_isolated import IsolatedStep from posydon.utils.posydonerror import FlowError +import copy import warnings @@ -88,24 +89,28 @@ def __call__(self,binary): merged_star_properties = self.merged_star_properties if self.verbose: print("Before Merger", binary.star_1.state,binary.star_2.state,binary.state, binary.event) - print("M1 , M2, he_core_mass1, he_core_mass2: ", binary.star_1.mass,binary.star_2.mass, binary.star_1.he_core_mass, binary.star_2.he_core_mass) - print("star_1.center_he4, star_2.center_he4, star_1.surface_he4, star_2.surface_he4: ", binary.star_1.center_he4,binary.star_2.center_he4, binary.star_1.surface_he4,binary.star_2.surface_he4) + print("M1 , M2, he_core_mass1, he_core_mass2, co_core_mass1, co_core_mass2: ", binary.star_1.mass,binary.star_2.mass, binary.star_1.he_core_mass, binary.star_2.he_core_mass, binary.star_1.co_core_mass, binary.star_2.co_core_mass) + print("star_1.center_h1, star_2.center_h1, star_1.center_he4, star_2.center_he4, star_1.center_c12, star_2.center_c12, star_1.surface_he4, star_2.surface_he4: ", binary.star_1.center_h1,binary.star_2.center_h1, binary.star_1.center_he4,binary.star_2.center_he4, binary.star_1.center_c12,binary.star_2.center_c12, binary.star_1.surface_he4,binary.star_2.surface_he4) if binary.state == "merged": if binary.event == 'oMerging1': binary.star_1,binary.star_2 = merged_star_properties(binary.star_1,binary.star_2) elif binary.event == 'oMerging2': binary.star_2,binary.star_1 = merged_star_properties(binary.star_2,binary.star_1) else: - raise FlowError("binary.state='merged' but binary.event != 'oMerging1/2'") + raise ValueError("binary.state='merged' but binary.event != 'oMerging1/2'") else: - raise FlowError("step_merging initiated but binary.state != 'merged'") + raise ValueError("step_merging initiated but binary.state != 'merged'") - binary.event = None if self.verbose: print("After Merger", binary.star_1.state,binary.star_2.state,binary.state, binary.event) - print("M_merged , he_core_mass merged: ", binary.star_1.mass, binary.star_1.he_core_mass) - print("star_1.center_he4, star_1.surface_he4: ", binary.star_1.center_he4, binary.star_1.surface_he4) + if binary.event == 'oMerging1': + print("M_merged , he_core_mass, co_core_mass merged: ", binary.star_1.mass, binary.star_1.he_core_mass, binary.star_1.co_core_mass) + print("center_h1, center_he4, center_c12, surface_h1, surface_he4: ", binary.star_1.center_h1, binary.star_1.center_he4, binary.star_1.center_c12, binary.star_1.surface_h1, binary.star_1.surface_he4) + elif binary.event == 'oMerging2': + print("M_merged , he_core_mass, co_core_mass merged: ", binary.star_2.mass, binary.star_2.he_core_mass, binary.star_2.co_core_mass) + print("center_h1, center_he4, center_c12, surface_h1, surface_he4: ", binary.star_2.center_h1, binary.star_2.center_he4, binary.star_2.center_c12, binary.star_2.surface_h1, binary.star_2.surface_he4) + binary.event = None super().__call__(binary) def merged_star_properties(self,star_base,comp): @@ -120,8 +125,8 @@ def merged_star_properties(self,star_base,comp): is the star that is engulfed """ #by default the stellar attributes that keep the same value from the - #merged_star = copy.copy(star_base) - merged_star = star_base + merged_star = copy.deepcopy(star_base) + #merged_star = star_base s1 = star_base.state s2 = comp.state @@ -143,6 +148,12 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma M2 = getattr(star2, "he_core_mass") - getattr(star2, "co_core_mass") else: M2 = getattr(star2, mass_weight2) + + if self.verbose: + print(abundance_name, mass_weight1,mass_weight2) + print("A_base, M_base_abundance, A_comp, M_comp_abundance", A1, M1, A2, M2) + print("mass_weighted_avg = ", (A1*M1 + A2*M2 ) / (M1+M2)) + return (A1*M1 + A2*M2 ) / (M1+M2) # MS + MS @@ -194,8 +205,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.log_LHe = star_base.log_LHe @@ -208,7 +218,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_HMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - merged_star = comp + merged_star = copy.deepcopy(comp) merged_star.mass = star_base.mass + comp.mass merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1", mass_weight2="H-rich_envelope_mass", mass_weight1="mass") @@ -222,8 +232,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.log_LHe = comp.log_LHe @@ -246,22 +255,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") #TODO : maybe he_core_mass makes more sense? merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: warnings.warn("weird compbination of CO core masses during merging") @@ -278,8 +293,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -298,8 +312,10 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has just a He core (is a HeMS star) pass # the central abundances are kept as the ones of star_base else: @@ -310,8 +326,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -320,7 +335,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_HeMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - merged_star = comp + merged_star = copy.deepcopy(comp) # add total and core masses for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: @@ -332,8 +347,10 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass == 0 and comp.co_core_mass > 0): # star_base is the HeMS Star and comp has a CO core pass # the central abundances are kept as the ones of star_base else: @@ -344,8 +361,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(star_base) @@ -364,22 +380,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: warnings.warn("weird compbination of CO core masses during merging") @@ -388,8 +410,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -398,7 +419,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_POSTHeMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - merged_star = comp + merged_star = copy.deepcopy(comp) # add total and core masses for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: current = getattr(merged_star, key) + getattr(star_base, key) @@ -408,22 +429,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: warnings.warn("weird compbination of CO core masses during merging") @@ -432,8 +459,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(star_base) @@ -452,22 +478,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: warnings.warn("weird compbination of CO core masses during merging") @@ -483,8 +515,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -496,21 +527,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma #WD is considered a stripped CO core # add total and core masses for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(merged_star, key) + getattr(comp, "mass") + if key not in ["c_core_mass", "o_core_mass"]: + current = getattr(merged_star, key) + getattr(comp, key) + else: + current = getattr(merged_star, key) + getattr(comp, "co_core_mass") setattr(merged_star, key,current) # weighted central abundances if merging cores. Else only from star_base - if (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has not + if (comp.co_core_mass is not None and star_base.co_core_mass == 0): # comp with CO core and the star_base has not merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 - elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): + merged_star.center_gamma = comp.center_gamma + elif (comp.co_core_mass is not None and star_base.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: warnings.warn("weird compbination of CO core masses during merging") @@ -519,8 +557,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -528,7 +565,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma # Star + NS/BH elif (s1 in STAR_STATES_NOT_CO and s2 in ["NS", "BH"]): - merged_star = comp + merged_star = copy.deepcopy(comp) # TODO: potentially flag a Thorne-Zytkov object massless_remnant = convert_star_to_massless_remnant(star_base) diff --git a/posydon/binary_evol/MESA/step_mesa.py b/posydon/binary_evol/MESA/step_mesa.py index 668aca0384..5938e5e020 100644 --- a/posydon/binary_evol/MESA/step_mesa.py +++ b/posydon/binary_evol/MESA/step_mesa.py @@ -252,7 +252,7 @@ def load_Interp(self, filename): """Load the interpolator that has been trained on the grid.""" if self.verbose: print("loading Interp: {}".format(filename)) - + # Check if interpolation files exist if not os.path.exists(filename): data_download() @@ -335,7 +335,7 @@ def __call__(self, binary): self.step(binary, interp_method=self.interpolation_method) if (self.stop_method == 'stop_at_max_time' and binary.time >= binary.properties.max_simulation_time): - + # self.flush_history = True # needed??? # stop_at_condition looks through the MESA output appended to the @@ -441,10 +441,10 @@ def update_properties_NN(self, star_1_CO=False, star_2_CO=False, setattr(binary, 'state', binary_state) setattr(binary, 'event', binary_event) setattr(binary, 'mass_transfer_case', MT_case) - + if binary.state == 'initial_RLOF': return - + if track_interpolation or self.save_initial_conditions: len_binary_hist = len(getattr(binary, "time_history")) length_binary_hist = len(cb_bh['age']) - 1 @@ -641,7 +641,7 @@ def update_properties_NN(self, star_1_CO=False, star_2_CO=False, getattr(star, key_h).extend(cb_bh[key_p][:-1]) elif key in ['state', 'lg_mdot']: continue - elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16']: + elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16', 'avg_c_in_c_core']: continue else: setattr(star, key, None) @@ -672,7 +672,7 @@ def update_properties_NN(self, star_1_CO=False, star_2_CO=False, setattr(binary, 'state', binary_state) setattr(binary, 'event', binary_event) setattr(binary, 'mass_transfer_case', MT_case) - + culmulative_mt_case = self.termination_flags[1] setattr(self.binary, f'culmulative_mt_case_{self.grid_type}', culmulative_mt_case) setattr(self.binary, f'interp_class_{self.grid_type}', interpolation_class) @@ -760,7 +760,7 @@ def update_properties_NN(self, star_1_CO=False, star_2_CO=False, history_of_attribute = (np.log10( 10**cb_bh[key_bh][0] + cf.bondi_hoyle( binary, accretor, donor, idx=len_binary_hist, - wind_disk_criteria=True, scheme='Kudritzki+2000'))) + wind_disk_criteria=True, scheme='Kudritzki+2000'))) if 10**history_of_attribute > edd: history_of_attribute = np.log10(edd) accretor.lg_mdot_history.append(history_of_attribute) @@ -886,7 +886,7 @@ def initial_final_interpolation(self, star_1_CO=False, star_2_CO=False): setattr(star, key, fv[key_p]) elif key == 'state': continue - elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16']: + elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16', 'avg_c_in_c_core']: continue else: setattr(star, key, None) @@ -896,7 +896,7 @@ def initial_final_interpolation(self, star_1_CO=False, star_2_CO=False): setattr(self.binary, f'interp_class_{self.grid_type}', interpolation_class) mt_history = self.classes['mt_history'] # mass transfer history (TF12 plot label) setattr(self.binary, f'mt_history_{self.grid_type}', mt_history) - + #TODO: add classifier for tf2 #setattr(self.binary, f'culmulative_mt_case', self.classes['termination_flags_2']) S1_state_inferred = cf.check_state_of_star(self.binary.star_1, @@ -933,7 +933,7 @@ def initial_final_interpolation(self, star_1_CO=False, star_2_CO=False): 10**fv[key_bh] + cf.bondi_hoyle( binary, accretor, donor, idx=-1, wind_disk_criteria=True, scheme='Kudritzki+2000')) - + mdot_edd = cf.eddington_limit(binary, idx=-1)[0] if 10**tmp_lg_mdot > mdot_edd: tmp_lg_mdot = np.log10(mdot_edd) @@ -1019,7 +1019,7 @@ def stop_at_condition(self, raise ValueError( 'Star can only be "star_1" or "star_2", you passed {0}'. format(star)) - + i = np.where(np.array(property_history) <= value)[0][-1] elif property in BINARYPROPERTIES: @@ -1067,7 +1067,7 @@ def stop_at_condition(self, continue interpolated_quanties[star][key] = self.interpolate_at_t( t, t_before, t_after, v_before, v_after) - + for key in BINARYPROPERTIES: if key in ['state', 'event', 'mass_transfer_case']: interpolated_quanties['binary'][key] = getattr( @@ -1258,7 +1258,7 @@ def __call__(self, binary): self.p_min <= p <= self.p_max): self.flip_stars_before_step = True super().__call__(self.binary) - + # redirect if outside grid for period elif (state_1 == 'H-rich_Core_H_burning' and state_2 == 'H-rich_Core_H_burning' and @@ -1394,7 +1394,7 @@ def __call__(self, binary): ecc == 0.)): super().__call__(self.binary) - + # period inside the grid, but m1 outside the grid elif ((not self.flip_stars_before_step and self.p_min <= p <= self.p_max and @@ -1498,7 +1498,7 @@ def __call__(self, binary): 'The star_1.state = %s, star_2.state = %s, binary.state = %s, ' 'binary.event = %s and not CO - HeMS - oRLO1/oRLO2!' % (state_1, state_2, state, event)) - + # redirect if outside grids if ((not self.flip_stars_before_step and self.m1_min <= m1 <= self.m1_max and @@ -1510,7 +1510,7 @@ def __call__(self, binary): self.p_min <= p <= self.p_max and ecc == 0.)): super().__call__(self.binary) - + # period inside the grid, but m1 outside the grid elif ((not self.flip_stars_before_step and self.p_min <= p <= self.p_max and @@ -1519,7 +1519,7 @@ def __call__(self, binary): set_binary_to_failed(self.binary) raise GridError(f'The mass of m1 ({m1}) is outside the grid,' 'while the period is inside the grid.') - + # period inside the grid, but m2 outside the grid elif ((not self.flip_stars_before_step and self.p_min <= p <= self.p_max and From 60490f97e2865655422432e6767a4f785726b38e Mon Sep 17 00:00:00 2001 From: ezapartas Date: Fri, 24 May 2024 16:48:41 +0300 Subject: [PATCH 2/4] added a try/except case for the mass_weighted_avg when an abundance (or avg_c_in_c_core) is np.nan to return np.nan --- posydon/binary_evol/DT/step_merged.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/posydon/binary_evol/DT/step_merged.py b/posydon/binary_evol/DT/step_merged.py index ea4336023f..ebad31aed3 100644 --- a/posydon/binary_evol/DT/step_merged.py +++ b/posydon/binary_evol/DT/step_merged.py @@ -149,12 +149,18 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma else: M2 = getattr(star2, mass_weight2) + try: + mass_weighted_avg_value = (A1*M1 + A2*M2 ) / (M1+M2) + except TypeError: + mass_weighted_avg_value = np.nan + if self.verbose: print(abundance_name, mass_weight1,mass_weight2) print("A_base, M_base_abundance, A_comp, M_comp_abundance", A1, M1, A2, M2) - print("mass_weighted_avg = ", (A1*M1 + A2*M2 ) / (M1+M2)) + print("mass_weighted_avg = ", mass_weighted_avg_value) + - return (A1*M1 + A2*M2 ) / (M1+M2) + return mass_weighted_avg_value # MS + MS if ( s1 in LIST_ACCEPTABLE_STATES_FOR_HMS From 5325335ff4309c7b45149b8ba6fdade81c5f9108 Mon Sep 17 00:00:00 2001 From: ezapartas Date: Thu, 6 Jun 2024 17:34:12 +0300 Subject: [PATCH 3/4] Update step_merged.py --- posydon/binary_evol/DT/step_merged.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/posydon/binary_evol/DT/step_merged.py b/posydon/binary_evol/DT/step_merged.py index ebad31aed3..e305262bd5 100644 --- a/posydon/binary_evol/DT/step_merged.py +++ b/posydon/binary_evol/DT/step_merged.py @@ -97,9 +97,9 @@ def __call__(self,binary): elif binary.event == 'oMerging2': binary.star_2,binary.star_1 = merged_star_properties(binary.star_2,binary.star_1) else: - raise ValueError("binary.state='merged' but binary.event != 'oMerging1/2'") + raise FlowError("binary.state='merged' but binary.event != 'oMerging1/2'") else: - raise ValueError("step_merging initiated but binary.state != 'merged'") + raise FlowError("step_merging initiated but binary.state != 'merged'") if self.verbose: print("After Merger", binary.star_1.state,binary.star_2.state,binary.state, binary.event) From 93ca5aa4036ca6a18249f579a4352a94d00ff82c Mon Sep 17 00:00:00 2001 From: ezapartas Date: Thu, 6 Jun 2024 17:40:44 +0300 Subject: [PATCH 4/4] Update step_merged.py --- posydon/binary_evol/DT/step_merged.py | 1 - 1 file changed, 1 deletion(-) diff --git a/posydon/binary_evol/DT/step_merged.py b/posydon/binary_evol/DT/step_merged.py index e305262bd5..0604dcf780 100644 --- a/posydon/binary_evol/DT/step_merged.py +++ b/posydon/binary_evol/DT/step_merged.py @@ -126,7 +126,6 @@ def merged_star_properties(self,star_base,comp): """ #by default the stellar attributes that keep the same value from the merged_star = copy.deepcopy(star_base) - #merged_star = star_base s1 = star_base.state s2 = comp.state