Skip to content

Commit

Permalink
Merge pull request #1280 from TeamCOMPAS/ge_qcrit_table_He
Browse files Browse the repository at this point in the history 
Ge qcrit table He
  • Loading branch information
@jeffriley
jeffriley authored Nov 22, 2024
2 parents 160e62d + d920139 commit 88f5a88
Show file tree
Hide file tree
Showing 14 changed files with 375 additions and 172 deletions.
6 changes: 3 additions & 3 deletions online-docs/pages/User guide/Program options/program-options-list-defaults.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,11 +311,11 @@ Default = 1.28

**--critical-mass-ratio-prescription** |br|
Critical mass ratio stability prescription (if any). |br|
Options: { NONE, ZERO, CLAEYS, GE20, GE20_IC, HURLEY_HJELLMING_WEBBINK } |br|
Options: { NONE, ZERO, CLAEYS, GE, GE_IC, HURLEY_HJELLMING_WEBBINK } |br|
``NONE`` defaults to the zeta prescription for stability. |br|
``CLAEYS`` uses qCrit values from Claeys et al. 2014. |br|
``GE20`` uses qCrit values from Ge et al. 2020 (adiabatic assumption). |br|
``GE20_IC`` uses qCrit values from Ge et al. 2020 (isentropic envelope assumption). |br|
``GE`` uses qCrit values from Ge et al. series (Papers I-V) (adiabatic assumption). |br|
``GE_IC`` uses qCrit values from Ge et al. series (Papers I-V) (isentropic envelope assumption). |br|
``HURLEY_HJELLMING_WEBBINK`` uses qCrit values from Hurley et al. 2002 (Hjellming & Webbink 1987 for mass transfer from a giant primary). |br|
Warning: if running with ``--critical-mass-ratio-prescription``, zetas will not be computed, so should not be trusted in the outputs. |br|
Default = NONE |br|
Expand Down
7 changes: 7 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,13 @@ What's new

Following is a brief list of important updates to the COMPAS code. A complete record of changes can be found in the file ``changelog.h``.

**03.08.02 Nov 18, 2024**

Updated implementation of the mass transfer stability critical mass ratio tables from the team of Hongwei Ge.
We now use all the most up to do date tables they've produced (public or otherwise) from Papers I-V, including
variations for adiabatic and isentropic treatments, variable accretion efficiency, and two different metallicities,
as well as for He stars (albeit only for adiabatic, fully conservative, solar metallicity stars).

**03.08.00 Nov 18, 2024**

Improved the treatment of stellar rotation (with further corrections in 03.08.01):
Expand Down
152 changes: 3 additions & 149 deletions src/BaseStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1399,9 +1399,9 @@ double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, c
qCrit = 0.0;
break;

case QCRIT_PRESCRIPTION::GE20:
case QCRIT_PRESCRIPTION::GE20_IC:
qCrit = CalculateCriticalMassRatioGe20(OPTIONS->QCritPrescription(), p_massTransferEfficiencyBeta);
case QCRIT_PRESCRIPTION::GE:
case QCRIT_PRESCRIPTION::GE_IC:
qCrit = CalculateCriticalMassRatioGeEtAl(OPTIONS->QCritPrescription(), p_massTransferEfficiencyBeta);
break;

case QCRIT_PRESCRIPTION::CLAEYS:
Expand All @@ -1427,152 +1427,6 @@ double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, c
}


/*
* Interpolate Ge+20 Critical Mass Ratios
*
* Function takes input QCRIT_PRESCRIPTION, currently either of the prescriptions for critical mass ratios
* from Ge et al. (2020), GE20 or GE20_IC. The first is the full adiabatic response, the second assumes
* artificially isentropic envelopes. From private communication with Ge, we have an updated datatable that
* includes qCrit for fully conservative and fully non-conservative MT, so we now interpolate on those as well.
*
* Interpolation is done linearly in logM, logR, and logZ
*
* double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, const double p_massTransferEfficiencyBeta)
*
* @param [IN] p_qCritPrescription Adopted critical mass ratio prescription
* @param [IN] p_massTransferEfficiencyBeta Mass transfer accretion efficiency
* @return Interpolated value of either the critical mass ratio or zeta for given stellar mass / radius
*/
double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, const double p_massTransferEfficiencyBeta) {

// Iterate over the two QCRIT_GE tables to get the qcrits at each metallicity
double qcritPerMetallicity[2];
std::vector<GE_QCRIT_TABLE> qCritTables = { QCRIT_GE_LOW_Z, QCRIT_GE_HIGH_Z };

for (int ii=0; ii<2; ii++) { // iterate over the vector of tables to store qCrits per metallicity

// Get vector of masses from qCritTable
GE_QCRIT_TABLE &qCritTable = qCritTables[ii];
DBL_VECTOR massesFromQCritTable = std::get<0>(qCritTable);
GE_QCRIT_RADII_QCRIT_VECTOR radiiQCritsFromQCritTable = std::get<1>(qCritTable);

INT_VECTOR indices = utils::BinarySearch(massesFromQCritTable, m_Mass);
int lowerMassIndex = indices[0];
int upperMassIndex = indices[1];

if (lowerMassIndex == -1) { // if masses are out of range, set to endpoints
lowerMassIndex = 0;
upperMassIndex = 1;
}
else if (upperMassIndex == -1) {
lowerMassIndex = massesFromQCritTable.size() - 2;
upperMassIndex = massesFromQCritTable.size() - 1;
}

// Get vector of radii from qCritTable for the lower and upper mass indices
std::vector<double> logRadiusVectorLowerMass = std::get<0>(radiiQCritsFromQCritTable[lowerMassIndex]);
std::vector<double> logRadiusVectorUpperMass = std::get<0>(radiiQCritsFromQCritTable[upperMassIndex]);

// Get the qCrit vector for the lower and upper mass bounds
std::vector<double> qCritVectorUpperEffLowerMass;
std::vector<double> qCritVectorUpperEffUpperMass;
std::vector<double> qCritVectorLowerEffLowerMass;
std::vector<double> qCritVectorLowerEffUpperMass;

// Set the appropriate qCrit vector, depends on MT eff and whether you use GE20 STD or IC
if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE20) {
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]);

}
}
else if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE20_IC) {
if (p_massTransferEfficiencyBeta > 0.5) {
qCritVectorUpperEffLowerMass = std::get<4>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorUpperEffUpperMass = std::get<4>(radiiQCritsFromQCritTable[upperMassIndex]);
qCritVectorLowerEffLowerMass = std::get<5>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorLowerEffUpperMass = std::get<5>(radiiQCritsFromQCritTable[upperMassIndex]);
}
else {
qCritVectorUpperEffLowerMass = std::get<5>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorUpperEffUpperMass = std::get<5>(radiiQCritsFromQCritTable[upperMassIndex]);
qCritVectorLowerEffLowerMass = std::get<6>(radiiQCritsFromQCritTable[lowerMassIndex]);
qCritVectorLowerEffUpperMass = std::get<6>(radiiQCritsFromQCritTable[upperMassIndex]);

}
}

// Get vector of radii from qCritTable for both lower and upper masses
INT_VECTOR indicesR0 = utils::BinarySearch(logRadiusVectorLowerMass, log10(m_Radius));
int lowerRadiusLowerMassIndex = indicesR0[0];
int upperRadiusLowerMassIndex = indicesR0[1];

if (lowerRadiusLowerMassIndex == -1) { // if radii are out of range, set to endpoints
lowerRadiusLowerMassIndex = 0;
upperRadiusLowerMassIndex = 1;
}
else if (upperRadiusLowerMassIndex == -1) {
lowerRadiusLowerMassIndex = logRadiusVectorLowerMass.size() - 2;
upperRadiusLowerMassIndex = logRadiusVectorLowerMass.size() - 1;
}

INT_VECTOR indicesR1 = utils::BinarySearch(logRadiusVectorUpperMass, log10(m_Radius));
int lowerRadiusUpperMassIndex = indicesR1[0];
int upperRadiusUpperMassIndex = indicesR1[1];

if (lowerRadiusUpperMassIndex == -1) { // if radii are out of range, set to endpoints
lowerRadiusUpperMassIndex = 0;
upperRadiusUpperMassIndex = 1;
}
else if (upperRadiusUpperMassIndex == -1) {
lowerRadiusUpperMassIndex = logRadiusVectorUpperMass.size() - 2;
upperRadiusUpperMassIndex = logRadiusVectorUpperMass.size() - 1;
}

// Set the 4 boundary points for the 2D interpolation
double qUppLowLow = qCritVectorUpperEffLowerMass[lowerRadiusLowerMassIndex];
double qUppLowUpp = qCritVectorUpperEffLowerMass[upperRadiusLowerMassIndex];
double qUppUppLow = qCritVectorUpperEffUpperMass[lowerRadiusUpperMassIndex];
double qUppUppUpp = qCritVectorUpperEffUpperMass[upperRadiusUpperMassIndex];
double qLowLowLow = qCritVectorLowerEffLowerMass[lowerRadiusLowerMassIndex];
double qLowLowUpp = qCritVectorLowerEffLowerMass[upperRadiusLowerMassIndex];
double qLowUppLow = qCritVectorLowerEffUpperMass[lowerRadiusUpperMassIndex];
double qLowUppUpp = qCritVectorLowerEffUpperMass[upperRadiusUpperMassIndex];

double logLowerMass = log10(massesFromQCritTable[lowerMassIndex]);
double logUpperMass = log10(massesFromQCritTable[upperMassIndex]);

double lowerLogRadiusLowerMass = logRadiusVectorLowerMass[lowerRadiusLowerMassIndex];
double upperLogRadiusLowerMass = logRadiusVectorLowerMass[upperRadiusLowerMassIndex];
double lowerLogRadiusUpperMass = logRadiusVectorUpperMass[lowerRadiusUpperMassIndex];
double upperLogRadiusUpperMass = logRadiusVectorUpperMass[upperRadiusUpperMassIndex];

// Interpolate on logR first, then logM, then on the mass transfer efficiency beta
double qCritUpperEffLowerMass = qUppLowLow + (upperLogRadiusLowerMass - log10(m_Radius)) / (upperLogRadiusLowerMass - lowerLogRadiusLowerMass) * (qUppLowUpp - qUppLowLow);
double qCritUpperEffUpperMass = qUppUppLow + (upperLogRadiusUpperMass - log10(m_Radius)) / (upperLogRadiusUpperMass - lowerLogRadiusUpperMass) * (qUppUppUpp - qUppUppLow);
double qCritLowerEffLowerMass = qLowLowLow + (upperLogRadiusLowerMass - log10(m_Radius)) / (upperLogRadiusLowerMass - lowerLogRadiusLowerMass) * (qLowLowUpp - qLowLowLow);
double qCritLowerEffUpperMass = qLowUppLow + (upperLogRadiusUpperMass - log10(m_Radius)) / (upperLogRadiusUpperMass - lowerLogRadiusUpperMass) * (qLowUppUpp - qLowUppLow);

double interpolatedQCritUpperEff = qCritUpperEffLowerMass + (logUpperMass - log10(m_Mass)) / (logUpperMass - logLowerMass) * (qCritUpperEffUpperMass - qCritUpperEffLowerMass);
double interpolatedQCritLowerEff = qCritLowerEffLowerMass + (logUpperMass - log10(m_Mass)) / (logUpperMass - logLowerMass) * (qCritLowerEffUpperMass - qCritLowerEffLowerMass);

double interpolatedQCritForZ = p_massTransferEfficiencyBeta * interpolatedQCritUpperEff + (1.0 - p_massTransferEfficiencyBeta) * interpolatedQCritLowerEff;
qcritPerMetallicity[ii] = interpolatedQCritForZ;
}
double logZlo = -3; // log10(0.001)
double logZhi = LOG10_ZSOL; // log10(0.02)

return qcritPerMetallicity[1] + (m_Log10Metallicity - logZhi)*(qcritPerMetallicity[1] - qcritPerMetallicity[0])/(logZhi - logZlo);
}


/*
Expand Down
7 changes: 4 additions & 3 deletions src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,8 +227,8 @@ class BaseStar {
virtual double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate,
const double p_massTransferEfficiencyBeta);
virtual double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const { return 0.0; } // Default is 0.0
double CalculateCriticalMassRatioGe20(const QCRIT_PRESCRIPTION p_qCritPrescription,
const double p_massTransferEfficiencyBeta) { return InterpolateGe20QCrit(p_qCritPrescription, p_massTransferEfficiencyBeta); }
virtual double CalculateCriticalMassRatioGeEtAl(const QCRIT_PRESCRIPTION p_qCritPrescription,
const double p_massTransferEfficiencyBeta) { return InterpolateGeEtAlQCrit(p_qCritPrescription, p_massTransferEfficiencyBeta); }
virtual double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.0; } // Default is 0.0

double CalculateDynamicalTimescale() const { return CalculateDynamicalTimescale_Static(m_Mass, m_Radius); } // Use class member variables
Expand Down Expand Up @@ -312,7 +312,8 @@ class BaseStar {

void HaltWinds() { m_Mdot = 0.0; } // Disable wind mass loss in current time step (e.g., if star is a donor or accretor in a RLOF episode)

double InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, const double p_massTransferEfficiencyBeta);
virtual double InterpolateGeEtAlQCrit(const QCRIT_PRESCRIPTION p_qCritPrescription,
const double p_massTransferEfficiencyBeta) { return 0.0; } // Placeholder, use interpolator for either H-rich or H-poor stars

void ResetEnvelopeExpulsationByPulsations() { m_EnvelopeJustExpelledByPulsations = false; }

Expand Down
Loading

0 comments on commit 88f5a88

Please sign in to comment.