Skip to content

[fix] Corrected/enhanced merger process#785

Open
ezapartas wants to merge 33 commits intomainfrom
manos_merger_copiedattributes
Open

[fix] Corrected/enhanced merger process#785
ezapartas wants to merge 33 commits intomainfrom
manos_merger_copiedattributes

Conversation

@ezapartas
Copy link
Contributor

Implementing and solving issues already discussed in a previous PR version #311 which is now deprecated.

  • No deepcopy during merger (which was working, but led to memory and speed issues)
  • However, correct caclulations of the mass weighted averages of abundances BEFORE mass changes.
  • calculating center_gamma and avg_c_in_c_core for mergers whenever possible (as the core value if the cores do not merge and np.nap and mass-weighted average respectively, if the cores merge).
  • added co_core_mass as a potential option in the matching criteria at track_match (mostly for mergers)

Reminder : The resulting merger product is continuing as a singlestar instance at the place of the engulfing star (i.e. if oMerging1, then star_1 would be the merged product) UNLESS a NS/BH is involved, where the latter is the only one surviving (potential Thorne-Zytkov object).

@ezapartas ezapartas requested review from maxbriel and sgossage and removed request for sgossage January 15, 2026 14:39
try:
mass_weighted_avg_value=(A1*M1+A2*M2)/(M1+M2)
except TypeError:
mass_weighted_avg_value= np.nan
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is setting the mass weighted average to np.nan going to cause issues later down the pipeline?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it possible to avoid a try/except?
Is it instead possible to verify the type of the values used?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  1. @astro-abhishek the np.nan may cause issues later, but I think is the most pythonic way of saying "the function was called byt physically there are issues at the calculation so we have no computed result, we don't know the value". If the question is whether it will make it crash, I am not sure, I think not, as np.nan is accepted as a value in calculations, it will just produce more "unknown" np.nan quantities wherever it is used in a calcualation.

  2. @maxbriel The following is a different way more exact checking, instead of try-except loop. Is this what you are thinking? If did not add it yet, but if you agree, please just add it as a commit directly.

den = M1 + M2
if not (np.isfinite(A1) and np.isfinite(A2) and np.isfinite(M1) and np.isfinite(M2)):
mass_weighted_avg_value = np.nan
elif den == 0:
mass_weighted_avg_value = np.nan
else:
mass_weighted_avg_value = (A1M1 + A2M2) / den


In both cases, if you have better coding ideas for this part (or any other) feel free to add them directly.

"log_R",
"center_c12"
"center_c12",
"co_core_mass"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what if the co_core_mass doesn't exist? (there's a comment about something similar right above this)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, this list assumes that the single star models and the respective EEPs have this column. If a user uses their one sinlge star and EEPs, not run with POSYDON, so not calculating a co_core (having only a c_core or an o_core for example), there may be a mismatch. For our POSYDON runs this is ok for now. If we want to generalize it so much, I think much of the structure should change anyway. I would leave it for now, keeping the above comment there for future reference.

@maxbriel maxbriel changed the title Corrected/enhanced merger process [fix] Corrected/enhanced merger process Feb 9, 2026
@maxbriel
Copy link
Collaborator

maxbriel commented Feb 9, 2026

There are two testing binaries in the suite that change with this PR.
@ezapartas can you confirm that these changes are expected from the change in the PR?

image image


# 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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is there this change from values to None? It happens later too.

Are there specific binaries that crash because the co_core_mass is undefined when reaching this stage?

Same n line 575.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, I don' remember exactly the reasoning. I think I thought that instead of raising a TypeError in case of comp.co_core_mass = None and crashing, I preferred it to go to the else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning").

But again, if you think of a better handling it , feel free to directly change it back

#these stellar attributes change value
merged_star.mass = (star_base.mass + comp.mass) * (1.-self.rel_mass_lost_HMS_HMS)

#TODO for key in ["center_h1", "center_he4", "center_c12", "center_n14","center_o16"]:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this TODO still needed?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, you are right, I took the line out

for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]:
current = getattr(star_base, key) + getattr(comp, key)
setattr(merged_star, key,current)
merged_star = comp
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are we changing what the merged star is here?

This is different than the definition of which one starts the RLO. This will make the companion the remaining star.

Copy link
Contributor Author

@ezapartas ezapartas Feb 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, actually this occurs even at line 256. Indeed, it "violates" the engulfind star (that starts the RLO) is our base of the merger, making the base to be the star with the most layers. So for example in a postMS+HMS merger, the postMS star becomes the base always. Same between postMS+HeMS. Even if it is not the postMS that engulfs. This is because the final merged star with have a more complex layering (envelope, core, etc). Of course in most cases it will indeed be the postMS star that is also initiating the RLO and the engulfing, but we allow for the other case too, giving the same outcome.

If needed we can restructure the code to remake the calculation using the engulfing always as the base, but this is only for the sake of clarity. It will not change the merged outcome in the end.

summary: The outcome in the end should not depend on which star engulfs (unless omeone argues otherwise). it is just easier to start the merger construction from that star in most cases, as it usually is the more "complex" star in the perspective of layers, i.e. it is a giant star

@maxbriel maxbriel added enhancement New feature or request bug Something isn't working labels Feb 10, 2026
@ezapartas
Copy link
Contributor Author

@astro-abhishek @maxbriel Thanks for the detailed comments. I tried to answer all of them and address whichever where doable. Feel free to directly change the coding way if you think a better option exists.

Note as soon as this PR is accepted, it is better to be pulled by the #788 before that PR is evaluated and merged.

@maxbriel
Copy link
Collaborator

Here's the output of the branch currently:
manos merger copiedattributes

main

The diff: https://www.diffchecker.com/7IurVLUA/


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)
Copy link
Contributor

@sgossage sgossage Mar 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we say *_base or e.g., *_1 or *_pri? I think other parts of the code use one of the latter. Also, should M_*_abundance be something like M_*_mass because M1/M2 hold masses?

Comment on lines +279 to +294
# companion with the envelope of star_base
parameters_to_mix = ["surface_h1", "surface_he4", "surface_c12", "surface_n14", "surface_o16"]
for abundance_name in parameters_to_mix:
setattr(merged_star, abundance_name, self.mass_weighted_avg(star_base,
comp,
abundance_name=abundance_name,
mass_weight1="H-rich_envelope_mass",
mass_weight2="mass"))

# The change of masses occurs after the calculation of weighted averages
# Note that the he core mass is unchanged, which means that
# adding the mass changes the envelope mass only.
merged_star.mass = star_base.mass + comp.mass

for key in STARPROPERTIES:
# these stellar attributes become np.nan
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"]:
setattr(merged_star, key, np.nan)
# Set parameters that are not expected to be meaningful after a merger to np.nan
self._apply_nan_attributes(merged_star)
Copy link
Contributor

@sgossage sgossage Mar 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some of this is repeated code in this if chain. I suggest generalizing this and creating a helper function for clarity/brevity that is called here and elsewhere, as appropriate.

Copy link
Contributor

@sgossage sgossage Mar 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

working on this here: sg_merger_copiedattributes

Copy link
Contributor

@sgossage sgossage Mar 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line (binary.event = None) means that the print statements based on binary.event (e.g., if binary.event == 'oMerging1':) below never fire when self.verbose = True. Should this be set to None after the verbose print? Is this line needed?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I commented the binary.event = None line for now.

"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!
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comment says that this for sure needs to be tested -- has it been?

@sgossage
Copy link
Contributor

sgossage commented Mar 21, 2026

I was working on generalizing the logic and creating functions to apply these merging rules and noticed something.

For post-HMS mergers, we have this logic (in step_merged.py):

                      ...
385                # star_base with CO core and the comp has a He core
386                elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0):
387                    pass # the central abundances are kept as the ones of star_base
                      ...

whereas for HMS-HMS mergers, we do not apply the same rule. Should we?

The following binary evolution from our test suite contains an HMS-HMS merger where the base star has a non-zero CO core mass while the companion has no CO core.

star_1= SingleStar(**{'mass': 9.845907 , 'state': 'H-rich_Core_H_burning','metallicity':1,
            'natal_kick_array':  [0.0, 4.190728383757787, 1.1521129697118118, 5.015343794234789]})
star_2 = SingleStar(**{'mass': 9.611029, 'state': 'H-rich_Core_H_burning','metallicity':1,
                'natal_kick_array': [0.0, 0.0, 0.0, 0.0]})
binary = BinaryStar(star_1, star_2, **{'time': 0.0,  'state': 'detached',  'event': 'ZAMS',
                        'orbital_period':  3.820571,'eccentricity': 0.0}, properties = sim_prop)

(The branch that I'm working on the logic in is sg_merger_copiedattributes. After we sort this out, if the changes there look good and the logic looks correct, we could merge it in to this branch. I've tested it with the test suite, but noticed that this binary produces a different outcome compared to manos_merger_copiedattribues because of this issue and how I currently wrote the logic.)

@maxbriel
Copy link
Collaborator

I was working on generalizing the logic and creating functions to apply these merging rules and noticed something.

For post-HMS mergers, we have this logic (in step_merged.py):

                      ...
385                # star_base with CO core and the comp has a He core
386                elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0):
387                    pass # the central abundances are kept as the ones of star_base
                      ...

whereas for HMS-HMS mergers, we do not apply the same rule. Should we?

The following binary evolution from our test suite contains an HMS-HMS merger where the base star has a non-zero CO core mass while the companion has no CO core.

star_1= SingleStar(**{'mass': 9.845907 , 'state': 'H-rich_Core_H_burning','metallicity':1,
            'natal_kick_array':  [0.0, 4.190728383757787, 1.1521129697118118, 5.015343794234789]})
star_2 = SingleStar(**{'mass': 9.611029, 'state': 'H-rich_Core_H_burning','metallicity':1,
                'natal_kick_array': [0.0, 0.0, 0.0, 0.0]})
binary = BinaryStar(star_1, star_2, **{'time': 0.0,  'state': 'detached',  'event': 'ZAMS',
                        'orbital_period':  3.820571,'eccentricity': 0.0}, properties = sim_prop)

(The branch that I'm working on the logic in is sg_merger_copiedattributes. After we sort this out, if the changes there look good and the logic looks correct, we could merge it in to this branch. I've tested it with the test suite, but noticed that this binary produces a different outcome compared to manos_merger_copiedattribues because of this issue and how I currently wrote the logic.)

This is strange. A HMS-HMS merger should not have a CO core or He core defined in either component, so the same rule shouldn't be applied.
It's concerning that something with a CO core is being put into the HMS-HMS category. A few clarifying questions: Is the stellar state being updated correctly during the evolution when entering the step? Is it still burning central hydrogen? And how large is the "CO core"?

@sgossage
Copy link
Contributor

sgossage commented Mar 22, 2026

It's concerning that something with a CO core is being put into the HMS-HMS category. A few clarifying questions: Is the stellar state being updated correctly during the evolution when entering the step? Is it still burning central hydrogen? And how large is the "CO core"?

@maxbriel this is the debugging output prior to merging:

Before Merger:
_____________

M1 = 4.212749495765598
M2 = 10.95210897562295
he_core_mass1 = 0.6365034990860114
he_core_mass2 = 2.019788928314233
co_core_mass1 = 0.10058851260570133
co_core_mass2 = 0.0

star_1.center_h1 = 0.03464608832551165
star_2.center_h1 = 0.016979750448212434
star_1.center_he4 = 0.8990294555960494
star_2.center_he4 = 0.9691717446871981
star_1.center_c12 = 0.05079687439526389
star_2.center_c12 = 9.913536056802633e-05
star_1.surface_he4 = 0.3547279661505059
star_2.surface_he4 = 0.32219608838969765

Star 1 is classed as H-rich_Core_H-burning by our criteria (log_LH = 3.4, surface_h1 = 0.6). I was able to trace it back to IF interpolation giving the star the 0.1 Msol CO core in step_HMS_HMS. The masses and age from IF interpolation seem reasonable, and it is classed as unstable MT (reasonable for that region of the grid), but somehow it determines a non-zero CO core mass, which is not reflected by our unstable MT models in the grids.

In the grids, this system may lie near a boundary between unstable MT and unstable reverse MT. While the unstable MT models have no CO core, the unstable reverse MT models do. These should be in different interpolation classes though?

I noticed that it seems we only have an unstable_MT interpolation class, and no unstable_reverse_MT. Could this be related to the issue? (@philipp-rajah)

At q = 0.95, the closest neighbors to this system are in the red boxed region:

q0 95_region

and at q = 0.99:

q1_region

(The axes in the plot are showing you the linear values of Porb and mass, not log10 values).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants