TeamCOMPAS / COMPAS

COMPAS rapid binary population synthesis code
http://compas.science
MIT License
65 stars 66 forks source link

BSE_Supernovae/Mass_Total@CO(SN) is occasionally zero #1033

Closed jeffriley closed 7 months ago

jeffriley commented 10 months ago

Describe the bug BSE_Supernovae/Mass_Total@CO(SN) is occasionally zero

Label the issue urgency_moderate - This is a moderately urgent issue severity_moderate - This is a moderately severe bug

To Reproduce Using COMPAS v02.41.02, run

./compas -n 1 --metallicity 0.02  --random-seed 1702719075

Expected behavior BSE_Supernovae/Mass_Total@CO(SN) should not be zero

Screenshots If applicable, add screenshots to help explain your problem.

Versioning (please complete the following information):

Additional context Add any other context about the problem here.

jeffriley commented 10 months ago

I believe this defect was introduced in this commit:

https://github.com/TeamCOMPAS/COMPAS/commit/4ca3ad4c88801b7ffb7a071a102637449435d545

as part of the work to add Accretion Induced Collapse (AIC) of ONeWD as another type of SN in COMPAS v02.31.08

Prior to the commit noted above, ONeWD::ResolveSupernova() called GiantBranch::ResolveSupernova(), which sets the variables

    m_SupernovaDetails.totalMassAtCOFormation  = m_Mass;
    m_SupernovaDetails.HeCoreMassAtCOFormation = m_HeCoreMass;
    m_SupernovaDetails.COCoreMassAtCOFormation = m_COCoreMass;
    m_SupernovaDetails.coreMassAtCOFormation   = m_CoreMass;

The commit noted above changed ONeWD::ResolveSupernova() to call ONeWD::ResolveAIC(). At some stage (I haven't checked where yet), this was refined so that in the current release (v02.41.02), WhiteDwarfs::ResolveSupernova() (called by all WDs) calls EvolveToNextPhase(). ONeWD::EvolveToNextPhase() calls ONeWD::ResolveAIC().

The bottom line is that ONeWD::ResolveAIC(), so AIC SN events, doesn't set the ...AtCOFormation variables that are set in GiantBranch::ResolveSupernova() - they are left at the default value of 0.0. They need to be set.

I don't know what the right solution is - I don't know what else in GiantBranch::ResolveSupernova() is also being missed that shouldn't be (I haven't looked). I don't know if ONEWD::ResolveAIC() can call GiantBranch::ResolveSupernova(), or whether we need to rewrite one or both so they can function together to set the variables i n question. Not sure duplicating code in ResolveAIC() is the right answer, but if it can't call (a modified form of) GiantBranch::ResolveSupernova() then maybe that's what we need to do.

@reinhold-willcox @ilyamandel

jeffriley commented 10 months ago

This is probably a defect for SNIA also (WhiteDwarfs::ResolveSNIA() also does not set the ...AtCOFormation variables).

reinhold-willcox commented 10 months ago

I'm happy to own this if I introduced the bug, however I am unable to reproduce it from the system given above. I still produce an AIC, but the MassTotal@CO is not 0, it's 1.38.

Screenshot from 2023-12-18 22-18-20

jeffriley commented 10 months ago

@reinhold-willcox What version of COMPAS are you using? Where do you get your numbers?

If I run this:

./compas -n 1 --metallicity 0.02  --random-seed 1702719075 --logfile-type csv

using COMPAS v02.41.02, I get the attached BSE_Supernovae file showing 0.0 for the @CO values.

BSE_Supernovae.csv

jeffriley commented 10 months ago

Ah - you have Record_Type in your output - you're pulling the values out of the detailed output file. The detailed output file is written each timestep (as you know), but the values that should be in the SN file are squirreled away in GiantBranch::ResolveSupernova() - and the new AIC (and SNIA) code bypasses that function. As noted above, I haven't checked if anything else from that function is also missed (that shouldn't be) because we don't call it for AIC and SNIA.

reinhold-willcox commented 10 months ago

This is BSE_Supernovae - we store Record_Type there too (apparently - it's in your csv file as well). That said, you're right that I was on an outdated branch, and I can reproduce the 0.0's in the different @CO masses with the latest dev.

However, the problem does not exist on somewhat more recent commits (more recent than the one where AIC and SNIA were introduced). I can compile on the commit starting 5fb489bf, which Mike and I worked on and which was committed on Dec. 13th. The bug is introduced afterwards, so I'm worried it's related to all of the reshuffling around the winds.

jeffriley commented 10 months ago

This is BSE_Supernovae - we store Record_Type there too (apparently - it's in your csv file as well).

Oh yeah - all files have record type :-) We only use it in the detailed output file and the pulsars file. Bad assumption on my part...

The bug is introduced afterwards, so I'm worried it's related to all of the reshuffling around the winds.

Oh? Well that's a problem. As I noted above, at some stage the calling pattern changed. I'm pretty sure that m_SupernovaDetails.totalMassAtCOFormation is only ever set in GiantBranch::ResolveSupernova(), which I don't think is (currently) called by the AIC and SNIA code - and it wasn't when that code was first introduced. Was a call added later?

jeffriley commented 10 months ago

@reinhold-willcox at first I thought my reapplication of the winds merge must have caused this problem, but I just checked the dev branch at the original winds commit (https://github.com/TeamCOMPAS/COMPAS/commit/e39c6b71173d7c1ddf404c26bf9c87a4376c4ab3), and the ...@CO columns in the BSE SN file are 0 there - so the problem was introduced at or before that commit.

But your commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7) was after the winds commit - and the code at your commit does not exhibit this problem. Your commit must have had a different parent from the winds commit?

EDIT: I wonder if we should freeze dev until we figure out what's happened?

reinhold-willcox commented 10 months ago

Ya something is looking really off here. We should take this offline and try to figure out where these issues cropped up. I'm hesitant to suggest reverting on the dev branch, that could cause a lot of headaches going forward. Maybe as a short term fix, we could spin up a new branch off an older version of dev, something like last-stable-release or something, and send out an announcement to anyone planning to do big runs in the short term.

jeffriley commented 10 months ago

Yeah, that's probably a good idea - but first we have to find a point where we are sure everything is good, then work forward. I will probably have some time over Christmas to do that if you don't. Let me know.

jeffriley commented 10 months ago

@reinhold-willcox @ilyamandel I'm now not convinced there is actually a problem with the repo.

We assumed there was because the problem for which this issue was opened is present in the original winds commit (https://github.com/TeamCOMPAS/COMPAS/commit/e39c6b71173d7c1ddf404c26bf9c87a4376c4ab3), but didn't appear to be present in commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7) (which was after the winds commit).

However, I think the problem is present in commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7), and in the current dev branch, but the original problem description didn't capture that - explanation follows:

The problem description indicates that this problem can be reproduced by executing

./compas -n 1 --metallicity 0.02 --random-seed 1702719075

but, as @reinhold-willcox pointed out, that does not reproduce the problem with the code at commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7) - thus creating the confusion.

The code at commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7) does not include the winds code - that was lost by the tides implementation merge, and reinstated after commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7). That is germane because without the winds code, the binary created from the run

./compas -n 1 --metallicity 0.02 --random-seed 1702719075

evolves differently from a run using (e.g.) the current dev branch that does include the winds code.

What's happening is that the binary created using the random seed 1702719075 evolves a star that becomes a TPAGB star, and when GiantBranch::ResolveSupernova() is called from BaseStar::UpdateAttributesAndAgeOneTimestep(), the call to TPAGB::IsSupernova() returns TRUE. The check in TPAGB::IsSupernova() is

utils::Compare(m_COCoreMass, m_GBParams[static_cast<int>(GBP::McSN)]) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0

and, in the case of seed 1702719075 returns TRUE because

m_COCoreMass = 1.38054, McSN = 1.38, and m_Mass = 1.38233

Because TPAGB::IsSupernova() returns TRUE in GiantBranch::ResolveSupernova(), the rest of the function is executed, and the ...@CO variables are set, and we don't see zeroes in the BSE_Supernovae file.

So under the circumstances (i.e. commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7) does not include the winds code) the problem description here was incomplete - the example given doesn't actually reproduce the problem.

However, by using the random seed 1703404893, by executing

./compas -n 1 --metallicity 0.02 --random-seed 1703404893

using the code at commit (https://github.com/TeamCOMPAS/COMPAS/commit/5fb489b4629637deb82e7a9a667899f6fa98a9a7), the problem is reproduced, and shows that what is happening is that, as before, the binary created evolves a star that becomes a TPAGB star, but this time when GiantBranch::ResolveSupernova() is called from BaseStar::UpdateAttributesAndAgeOneTimestep(), the call to TPAGB::IsSupernova() returns FALSE. Again, the check in TPAGB::IsSupernova() is

utils::Compare(m_COCoreMass, m_GBParams[static_cast<int>(GBP::McSN)]) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0

but in the case of seed 1703404893 it returns FALSE because

m_COCoreMass = 1.37984, McSN = 1.38, m_Mass = 6.84345

and because TPAGB::IsSupernova() returns FALSE in GiantBranch::ResolveSupernova(), the rest of the function is not executed, and the ...@CO variables are not set, so we see zeroes in the BSE_Supernovae file. The BSE_Supernovae entry for this binary is written because the star eventually evolves to a ONeWD, and ONeWD::EvolveToNextPhase() callse ResolveAIC() , which in turn calls IsMassAboveEcsnThreshold() which returns TRUE because

m_Mass = 1.38027, MECS = 1.38

So we've lost the envelope but m_Mass is now a little bigger than m_COCoreMass was earlier, and we pass the check - so we flag the SN and write the record to then BSE_Supernovae file, but since GiantBranch::ResolveSupernova() wasn't called earlier we get zeros for the ...@CO variables.

reinhold-willcox commented 10 months ago

Good catch! I'm glad the issue wasn't with the merges.

It looks like the solution is then finding a way to set the @CO variables for AIC and SNIA. I can take a look at this in the next few weeks. And the current dev should be OK for everyone to use again.

ilyamandel commented 10 months ago

Thank you for investigating, @jeffriley !

Just to make sure I understand correctly:

Is this correct?

jeffriley commented 10 months ago

@ilyamandel Yes, I believe that is correct - and confusion in how the problem presented caused us to think there may have been a problem with the repo, but better understanding the problem resolved that.

@reinhold-willcox Yes, but (always a but...), I think rather than just setting the ...@CO variables, the answer lies in working out why Giantbranch::ResolveSupernva() isn't being called (beyond the superficial reason that the check wasn't satisfied) - it seems that we don't think it is a SN in the usual code path, but then we do when we decide we have an AIC. Is that the way it should be, or is there something else we should be checking so that Giantbranch::ResolveSupernva() is called and the variables correctly set there?

jeffriley commented 10 months ago

@reinhold-willcox I think the following code snippet from GiantBranch::ResolveSupernova() is also being missed for any SN event that doesn't call GiantBranch::ResolveSupernova():

        if (OPTIONS->EvolutionMode() == EVOLUTION_MODE::SSE) {                                      // only if SSE (BSE does its own SN printing)
        StashSupernovaDetails(stellarType);
    }

I think missing this code means that the supernova event is not record the in SSE_Supernovae file when COMPAS is run in SSE mode. (See the comments in GiantBranch::ResolveSupernova() for an explanation of why the call to StashSupernovaDetails() is required)