Repository navigation
Grid interpolation for Ge's mass transfer stability model #1484
Description
The interpolation of the q_crit grids for mass transfer stability with Ge's model is giving weird results (--critical-mass-ratio-prescription GE or GE_IC).
severity_moderate
urgency_moderate
The figure below shows post first stable mass transfer binaries with a mass ratio q~1, thermal timescale in red, nuclear timescale in green. The donors are HG or MS stars. The vertical lines show the donor masses with q_crit values in the grid and match a saw tooth-like upper limit of the red points tracing where mass transfer gets unstable.
Lines 1320-26 in MainSequence.cpp, in MainSequence::InterpolateGeEtAlQCrit:
double interpolatedQCritUpperEff = (logMass < logLowerMass) ? qCritUpperEffLowerMass : (logMass > logUpperMass) ? qCritUpperEffUpperMass : qCritUpperEffLowerMass + (logUpperMass - logMass) / (logUpperMass - logLowerMass) * (qCritUpperEffUpperMass - qCritUpperEffLowerMass); double interpolatedQCritLowerEff = (logMass < logLowerMass) ? qCritLowerEffLowerMass : (logMass > logUpperMass) ? qCritLowerEffUpperMass : qCritLowerEffLowerMass + (logUpperMass - logMass) / (logUpperMass - logLowerMass) * (qCritLowerEffUpperMass - qCritLowerEffLowerMass);
The two numerators (logUpperMass - logMass) should be (logMass - logLowerMass). As they currently are, logMass = logUpperMass would give qCritLowerMass, and logMass = logLowerMass would give qCritUpperMass, with an inverted scaling in logMass explaining the saw tooth aspect of the figure.
If this comment about the mass is correct, there are similar issues regarding the interpolation of the radius and beta above and below the code interpolating the mass in the same function.
Versioning
COMPAS v03.29.05 (gsl v2.8, boost v1.92.0, HDF5 v2.2.0)
Activity
ilyamandel commentedon Aug 21, 2026
Hi @Dylan-Hebrail ,
Thank you for finding this issue and carefully investigating it!
@reinhold-willcox , I believe this is your code -- would you be willing to take a look at it?
reinhold-willcox commentedon Aug 24, 2026
Ya woof, that's on me. Thanks for catching this @Dylan-Hebrail ! Agreed that this same mistake is made in all the interpolations in this function.
Do you have time to address this in the short term? If not, I can try to get to it, but it may take a few weeks.
Dylan-Hebrail commentedon Aug 24, 2026
Hi @reinhold-willcox and thank you for your answer! I'm happy to address it in a couple of days and run some testing.
When reading the function in more details, I also noticed something for the beta interpolation. If I am correct these lines assign lower and upper qcrit values for
p_massTransferEfficiencyBeta> 0.5 or <= 0.5, meaning the interpolation is done between beta = 0 and 0.5 or beta = 0.5 and 1.`if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE) {
if (p_massTransferEfficiencyBeta > 0.5) {
qCritVectorUpperEffLowerMass = std::get<1>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorUpperEffUpperMass = std::get<1>(radiiQCritsFromQCritTable[upperMassIndex]);
qCritVectorLowerEffLowerMass = std::get<2>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorLowerEffUpperMass = std::get<2>(radiiQCritsFromQCritTable[upperMassIndex]);
}
else {
qCritVectorUpperEffLowerMass = std::get<2>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorUpperEffUpperMass = std::get<2>(radiiQCritsFromQCritTable[upperMassIndex]);
qCritVectorLowerEffLowerMass = std::get<3>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorLowerEffUpperMass = std::get<3>(radiiQCritsFromQCritTable[upperMassIndex]);
(lines 1227-1239)
However, towards the end of the function, the interpolation is done with this line:
double interpolatedQCritForZ = p_massTransferEfficiencyBeta * interpolatedQCritUpperEff + (1.0 - p_massTransferEfficiencyBeta) * interpolatedQCritLowerEff;(line 1327)Since
p_massTransferEfficiencyBetais not updated, it is still the value of beta and not an interpolation factor between the boundaries. For example ifp_massTransferEfficiencyBetais 0.5, the condition will interpolate between the grid beta = 0 and beta = 0.5, and give the value halfway between these two: the value for beta = 0.25.If you agree with this, I can fix it by defining a local interpolation factor within the beta interval.
ilyamandel commentedon Sep 24, 2026
Closing -- fixed in #1486