Divide-by-zero in EAGB::CalculateTimescales() #1348
Description
Describe the bug
A clear and concise description of what the bug is.
Running COMPAS v03.15.00 with:
./compas -n 1 --random-seed 0
results in a divide-by-zero error in EAGB::CalculateTimescales(), with the following stack trace:
EAGB::CalculateTimescales(double, std::vector<double, std::allocator<double> >&)
TPAGB::CalculateTimescales(double, std::vector<double, std::allocator<double> >&)
Remnants::CalculateTimescales(double, std::vector<double, std::allocator<double> >&)
Remnants::CalculateTimescales()
BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(double, double, double)
BaseStar::UpdateAttributesAndAgeOneTimestep(double, double, double, bool, bool)
Star::UpdateAttributesAndAgeOneTimestep(double, double, double, bool, bool)
Star::AgeOneTimestep(double, bool)
BaseBinaryStar::EvolveOneTimestep(double)
BaseBinaryStar::Evolve()
BinaryStar::Evolve()
EvolveBinaryStars()
main
The line of code in EAGB::CalculateTimescales() at which the divide-by-zero occurs is:
timescales(tinf1_FAGB) = tBAGB + ((1.0 / (p1 * gbParams(AHe) * gbParams(D))) * PPOW((gbParams(D) / LBAGB), p1_p));
The problem occurs because the variable LBAGB is 0.0, having been calculated at this line in EAGB::CalculateTimescales():
double LBAGB = CalculateLuminosityAtBAGB(p_Mass);
because p_Mass is passed to EAGB::CalculateTimescales() as 0.0.
The short explanation is that at some stage during the evolution of the binary, star 1 becomes a neutron star, and in the course of evolution we call BaseStar::UpdateAttributesAndAgeOneTimestep() which calls CalculateTimescales() which for the NS calls Remnants::CalculateTimescales(), which in turn calls TPAGB:CalculateTimescales() with m_Mass0 as the mass parameter, which calls EAGB::CalculateTimescales() with p_mass = 0.0, leading to the divide-by-zero error.
m_Mass0 is initialised to 0.0 in NS::Initialise(). Should it be? Not sure why... I don't think this was a change in the recent pulsar PR (#1315). It's possible it has been there for a while and we have not noticed it.
@ilyamandel, might it be a side-effect of one of your recent changes (I think I remember you rearranged some of the mass changes and maybe added/moved a call to one of the "UpdateAttributes" functions)? Even if it is, maybe we shouldn't be initialising m_Mass0 to 0.0 in NS::Initialise().
- OS: Ubuntu 22.04.2 LTS
- COMPAS v03.15.00
EDIT: I set the severity to major and urgency high because this defect sets one of the timescales to who-knows-what (probably +inf), and that value is used at every timestep to calculate various stellar attributes - and who knows what side-effects that has...