TeamCOMPAS / COMPAS

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

A bug in supernovae #1149

Open ilyamandel opened 4 months ago

ilyamandel commented 4 months ago

We have a bug in SSE [and possibly BSE -- I only checked SSE for now].

Example command: ./COMPAS --mode SSE --initial-mass 6.20 (with COMPAS v02.49.04)

This ultimately triggers a call to ResolveElectronCaptureSupernova(), which should therefore leave behind a OXYGEN_NEON_WHITE_DWARF [checked by putting in print statements].

However, that call only happens via ResolveMassLoss() in Star::EvolveOneTimestep(). And ResolveMassLoss() doesn't actually execute the switch to this new type -- instead, it throws a warning because it doesn't expect the type to change. So at the end of ResolveMassLoss(), from which ResolveElectronCaptureSupernova() was called, the type is still TPAGB [checked by putting in print statements].

Then, after EvolveOneTimestep() is called from Star::Evolve(), we immediately call UpdateAttributes(0.0. 0.0, true). By now, however, the star has lost its remaining envelope because it went through a SN, as ResolveElectronCaptureSupernova() set the core mass and the mass to the same value. Yet it's still a star of the TPAGB type, as its type and class have not yet been updated.

So when ResolveSupernova() is called from UpdateAttributes, TPAGB::IsSupernova() evaluates to false, and a SN doesn't happen. As a result, we are left with a COWD, not an ONeWD.

Now, one can quibble whether this should have been a COWD or an ONeWD in reality -- but I think it's quite clear that the code is not doing what was intended, i.e., not setting the star to be an ONeWD as the initial evaluation of ResolveElectronCaptureSupernova() demands.

My preferred solution would be for a star to actually go supernova and change type accordingly when it successfully goes through a supernova function such as ResolveElectronCaptureSupernova(). A supernova is a more "extreme" event than mere mass loss, and so should preempt other behaviour. However, since the order of calls here is very convoluted, I am rather afraid of how many things I might break if I try to fix this -- so, @jeffriley , could you please take a look?

jeffriley commented 4 months ago

I don't really like the idea of switching-in-place with the code the way it is. Hurley's BSE effectively is a big switch statement for stellar type, and he will change stellar type in one of the cases, but continue executing the remainder of the code in that case - so effectively executing code that is for the previous stellar type... We kind-of do the same thing, in that we decide we want to switch to a new stellar type, but only do it at the timestep boundary - with the exception (now) that we switch the secondary to a MR for merger products.

I don't really like our current approach (we finish the timestep as though we haven't switched, but we know we should switch...), but I also don't like switching-in-place - because once we switch we are (or could be) executing a function for the previous stellar type, and then we might end up calling a function that will now be for the new stellar type - we could end up in some weird state.

What I would like to do is allow switch-in-place, but once the switch is done, end the timestep and start a new one - so we always execute a timestep as just one stellar type. The problem with that - at least I think it may be a problem - is that we calculate mass loss, mass transfer etc. at the start of the timestep, based on the expected duration of the timestep, and if we end the timestep early what does that mean? What if when we switch we have already applied mass loss and mass transfer - do we end the timestep at the switch and assume it was the expected duration? Maybe we do... That was one of my worries about "only doing one thing during a timestep". I can see the attraction, but I think we have to be careful about how we actually implement that. And if we do end a timestep "early", what might we not have done to set up the next timestep? Maybe nothing, but something we should check

I think I can change the code so that we switch-in-place (i.e. when we decide we should switch, we do, then end the timestep), but I'd prefer to do that after the error handling PR is merged - it will be easier to do with the infrastructure that PR sets up. And once we have a mechanism for sensibly ending a timestep "early", we can use it for any significant event - not just switching stellar type.

ilyamandel commented 4 months ago

@jeffriley -- as you already know, I agree with you, so just publicly confirming this here.

As you also know, I think of your comment about needing to check what we might have forgotten to do when dealing with the change in stellar type as one more argument in favour of doing one "important" thing per time step. We want to ease the burden of future development, and that means limiting the range of possible implications that every developer has to think through. :-)

My current thinking is that we should check for dramatic things like SNe at the start of a time step. If the conditions are satisfied, we should take a size 0 step to take care of these special events. If not, we take a normal time step. But I agree that this can be deferred to be part of a broader change to time stepping.