-inf (divide-by-zero) in BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023()

[jeffriley]

Running COMPAS with:

./compas -n1 --minimum-secondary-mass 5

results in -inf in BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023() in the following statment:

double logMdotUncorrected  = log10(p_Mdot);     // uncorrected mass-loss rate

because the parameter p_Mdot = 0.0

double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const {

    const double teffRef = 141.0E3;                                 // reference effective temperature in Kelvin
    const double teffMin = 100.0E3;                                 // minimum effective temperature in Kelvin to apply correction

    double teff                = m_Temperature * TSOL;              // get effective temperature in Kelvin
    double logMdotUncorrected  = log10(p_Mdot);                     // uncorrected mass-loss rate
    double logMdotCorrected    = 0.0;

    // Only apply to sufficiently hot stars
    if (utils::Compare(teff, teffMin) > 0) {
        logMdotCorrected = logMdotUncorrected - 6.0 * log10(teff / teffRef);
    }
    else {
        logMdotCorrected = logMdotUncorrected;
    }

    return PPOW(10.0, logMdotCorrected);
}

The naive solution is to add the following sanity check as the first line of the function:

if (p_Mdot <= 0.0) return 0.0;                                  // don't use utils::Compare() here

so we just return a mass loss rate of 0.0 if the rate passed in p_Mdot is 0.0

p_Mdot can be passed as 0.0 from HeMS::CalculateMassLossRateMerritt2024() when CalculateMassLossRateWolfRayetSanderVink2020() is called with p_Mu < 1.0 (as it is in HeMS::CalculateMassLossRateMerritt2024()) (this is the case that caused the failure that prompted opening of this issue).
There may be other ways for the problem to manifest

@jmerritt1, @SimonStevenson, can you confirm that the suggestion above is the best solution?

To Reproduce
Run COMPAS with:

./compas -n1 --minimum-secondary-mass 5

Expected behavior
Production code does not cause -inf

**Versioning

• OS: Ubuntu 22.04
• COMPAS v03.07.01

[SimonStevenson]

It is true that the Sander+ 2023 prescription gives zero mass loss for low luminosity helium stars. However I thought that in our default prescription we apply a minimum mass loss for these stars from Vink2017. So I'm not sure why Mdot is zero.

[EDIT: ah we do, on line 422 in HeMS.cpp. So I think this issue should not have any real impact. Definitely not a major bug -- a major bug should result in incorrect results for a sizeable fraction of systems.]

Your suggestion:
return a mass loss rate of 0.0 if the rate passed in p_Mdot is 0.0
Sounds good to me.

[jeffriley]

Well, it depends on your definition of a major defect - I think dividing by zero always indicates a problem that should be fixed. As I mention elsewhere, apart from doing something we shouldn't (dividing by zero), it messes up the fp-error handling.

[ilyamandel]

@jeffriley , your example didn't tell me which random seed to use, so I didn't actually reproduce it -- but the solution you proposed seems sensible and @SimonStevenson approved it, so I implemented it in #1291 . Closing.