Skip to content

Commit

Permalink
added functionality for Ge et al. critical mass ratios for non-conser…
Browse files Browse the repository at this point in the history 
…vative MT
  • Loading branch information
@reinhold-willcox
reinhold-willcox committed Mar 29, 2024
1 parent 83b4f3c commit 5e5b8da
Show file tree
Hide file tree
Showing 6 changed files with 55 additions and 36 deletions.
4 changes: 1 addition & 3 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1890,10 +1890,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
else if (OPTIONS->QCritPrescription() != QCRIT_PRESCRIPTION::NONE) { // Determine stability based on critical mass ratios

// NOTE: Critical mass ratio is defined as mAccretor/mDonor
double qCrit = m_Donor->CalculateCriticalMassRatio(m_Accretor->IsDegenerate());

double qCrit = m_Donor->CalculateCriticalMassRatio(m_Accretor->IsDegenerate(), m_FractionAccreted);
isUnstable = (m_Accretor->Mass() / m_Donor->Mass()) < qCrit;
m_FractionAccreted = 1.0; // Accretion is assumed fully conservative for qCrit calculations
}
else { // Determine stability based on zetas

Expand Down
68 changes: 42 additions & 26 deletions src/BaseStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1324,20 +1324,20 @@ double BaseStar::CalculateZetaAdiabaticSPH(const double p_CoreMass) const {
/*
* Calculate the critical mass ratio for unstable mass transfer
*
* double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate)
* double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, const double p_massTransferEfficiencyBeta)
*
* @param [IN] p_AccretorIsDegenerate Whether or not the accretor is a degenerate star
* @return Critical mass ratio
*/
double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) {
double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, const double p_massTransferEfficiencyBeta) {

double qCrit = 0.0;
QCRIT_PRESCRIPTION qCritPrescription = OPTIONS->QCritPrescription();

switch (qCritPrescription) {
case QCRIT_PRESCRIPTION::GE20:
case QCRIT_PRESCRIPTION::GE20_IC:
qCrit = CalculateCriticalMassRatioGe20(qCritPrescription);
qCrit = CalculateCriticalMassRatioGe20(qCritPrescription, p_massTransferEfficiencyBeta);
break;
case QCRIT_PRESCRIPTION::CLAEYS:
qCrit = CalculateCriticalMassRatioClaeys14(p_AccretorIsDegenerate);
Expand All @@ -1358,18 +1358,19 @@ double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) {
*
* 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.
* 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.
*
* double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription)
* double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, const double p_massTransferEfficiencyBeta)
*
* @return Interpolated value of either the critical mass ratio or zeta for given stellar mass / radius
*/
double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription) {
double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, const double p_massTransferEfficiencyBeta) {

// Get vector of masses from GE20_QCRIT_AND_ZETA
std::vector<double> massesFromGe20 = std::get<0>(GE20_QCRIT_AND_ZETA);
// Get vector of masses from GE20_QCRIT
std::vector<double> massesFromGe20 = std::get<0>(GE20_QCRIT);
std::vector< std::tuple<std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>>>
radiiQCritsZetasFromGe20 = std::get<1>(GE20_QCRIT_AND_ZETA);
radiiQCritsZetasFromGe20 = std::get<1>(GE20_QCRIT);

std::vector<int> ind = utils::binarySearch(massesFromGe20, m_Mass);
int lowerMassInd = ind[0];
Expand All @@ -1384,25 +1385,31 @@ double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescripti
upperMassInd = massesFromGe20.size() - 1;
}

// Get vector of radii from GE20_QCRIT_AND_ZETA for the lower and upper mass indices
// Get vector of radii from GE20_QCRIT for the lower and upper mass indices
std::vector<double> logRadiusVectorLowerMass = std::get<0>(radiiQCritsZetasFromGe20[lowerMassInd]);
std::vector<double> logRadiusVectorUpperMass = std::get<0>(radiiQCritsZetasFromGe20[upperMassInd]);

// Get the qCrit vector for the lower and upper mass bounds
std::vector<double> qCritVectorLowerMass;
std::vector<double> qCritVectorUpperMass;
std::vector<double> qCritFullVectorLowerMass;
std::vector<double> qCritFullVectorUpperMass;
std::vector<double> qCritNoncVectorLowerMass;
std::vector<double> qCritNoncVectorUpperMass;

// One of the following must be set
// One of the following must be set - Full and Nonc distinguish fully conservative and non-conservative MT
if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE20) {
qCritVectorLowerMass = std::get<1>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritVectorUpperMass = std::get<1>(radiiQCritsZetasFromGe20[upperMassInd]);
qCritFullVectorLowerMass = std::get<1>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritFullVectorUpperMass = std::get<1>(radiiQCritsZetasFromGe20[upperMassInd]);
qCritNoncVectorLowerMass = std::get<3>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritNoncVectorUpperMass = std::get<3>(radiiQCritsZetasFromGe20[upperMassInd]);
}
else if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE20_IC) {
qCritVectorLowerMass = std::get<2>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritVectorUpperMass = std::get<2>(radiiQCritsZetasFromGe20[upperMassInd]);
qCritFullVectorLowerMass = std::get<2>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritFullVectorUpperMass = std::get<2>(radiiQCritsZetasFromGe20[upperMassInd]);
qCritNoncVectorLowerMass = std::get<4>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritNoncVectorUpperMass = std::get<4>(radiiQCritsZetasFromGe20[upperMassInd]);
}

// Get vector of radii from GE20_QCRIT_AND_ZETA for both lower and upper masses
// Get vector of radii from GE20_QCRIT for both lower and upper masses
std::vector<int> indR0 = utils::binarySearch(logRadiusVectorLowerMass, log10(m_Radius));
double lowerRadiusLowerMassInd = indR0[0];
double upperRadiusLowerMassInd = indR0[1];
Expand Down Expand Up @@ -1430,10 +1437,14 @@ double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescripti
}

// Set the 4 boundary points for the 2D interpolation
double qLowLow = qCritVectorLowerMass[lowerRadiusLowerMassInd];
double qLowUpp = qCritVectorLowerMass[upperRadiusLowerMassInd];
double qUppLow = qCritVectorUpperMass[lowerRadiusUpperMassInd];
double qUppUpp = qCritVectorUpperMass[upperRadiusUpperMassInd];
double qLowLowFull = qCritFullVectorLowerMass[lowerRadiusLowerMassInd];
double qLowUppFull = qCritFullVectorLowerMass[upperRadiusLowerMassInd];
double qUppLowFull = qCritFullVectorUpperMass[lowerRadiusUpperMassInd];
double qUppUppFull = qCritFullVectorUpperMass[upperRadiusUpperMassInd];
double qLowLowNonc = qCritNoncVectorLowerMass[lowerRadiusLowerMassInd];
double qLowUppNonc = qCritNoncVectorLowerMass[upperRadiusLowerMassInd];
double qUppLowNonc = qCritNoncVectorUpperMass[lowerRadiusUpperMassInd];
double qUppUppNonc = qCritNoncVectorUpperMass[upperRadiusUpperMassInd];

double lowerMass = massesFromGe20[lowerMassInd];
double upperMass = massesFromGe20[upperMassInd];
Expand All @@ -1443,11 +1454,16 @@ double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescripti
double lowerRadiusUpperMass = PPOW(10.0, logRadiusVectorUpperMass[lowerRadiusUpperMassInd]);
double upperRadiusUpperMass = PPOW(10.0, logRadiusVectorUpperMass[upperRadiusUpperMassInd]);

// Interpolate on the radii first, then the masses
double qCritLowerMass = qLowLow + (upperRadiusLowerMass - m_Radius) / (upperRadiusLowerMass - lowerRadiusLowerMass) * (qLowUpp - qLowLow);
double qCritUpperMass = qUppLow + (upperRadiusUpperMass - m_Radius) / (upperRadiusUpperMass - lowerRadiusUpperMass) * (qUppUpp - qUppLow);
double interpolatedQCrit = qCritLowerMass + (upperMass - m_Mass) / (upperMass - lowerMass) * (qCritUpperMass - qCritLowerMass);
// Interpolate on the radii first, then the masses, then on the mass transfer efficiency beta
double qCritFullLowerMass = qLowLowFull + (upperRadiusLowerMass - m_Radius) / (upperRadiusLowerMass - lowerRadiusLowerMass) * (qLowUppFull - qLowLowFull);
double qCritFullUpperMass = qUppLowFull + (upperRadiusUpperMass - m_Radius) / (upperRadiusUpperMass - lowerRadiusUpperMass) * (qUppUppFull - qUppLowFull);
double qCritNoncLowerMass = qLowLowNonc + (upperRadiusLowerMass - m_Radius) / (upperRadiusLowerMass - lowerRadiusLowerMass) * (qLowUppNonc - qLowLowNonc);
double qCritNoncUpperMass = qUppLowNonc + (upperRadiusUpperMass - m_Radius) / (upperRadiusUpperMass - lowerRadiusUpperMass) * (qUppUppNonc - qUppLowNonc);

double interpolatedQCritFull = qCritFullLowerMass + (upperMass - m_Mass) / (upperMass - lowerMass) * (qCritFullUpperMass - qCritFullLowerMass);
double interpolatedQCritNonc = qCritNoncLowerMass + (upperMass - m_Mass) / (upperMass - lowerMass) * (qCritNoncUpperMass - qCritNoncLowerMass);

double interpolatedQCrit = p_massTransferEfficiencyBeta * interpolatedQCritNonc + (1-p_massTransferEfficiencyBeta)*interpolatedQCritFull;
return interpolatedQCrit;
}

Expand Down
8 changes: 5 additions & 3 deletions src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,11 @@ class BaseStar {

virtual double CalculateConvectiveEnvelopeMass() const { return 0.0; }

virtual double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate);
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) { return InterpolateGe20QCrit(p_qCritPrescription); }
double CalculateCriticalMassRatioGe20(const QCRIT_PRESCRIPTION p_qCritPrescription,
const double p_massTransferEfficiencyBeta) { return InterpolateGe20QCrit(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 @@ -228,7 +230,7 @@ class BaseStar {

virtual void ResolveShellChange(const double p_AccretedMass) { } // Default does nothing, use inheritance for WDs.

double InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription);
double InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, const double p_massTransferEfficiencyBeta);

void ResolveAccretion(const double p_AccretionMass) { m_Mass = std::max(0.0, m_Mass + p_AccretionMass); } // Handles donation and accretion - won't let mass go negative

Expand Down
3 changes: 2 additions & 1 deletion src/Remnants.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@ class Remnants: virtual public BaseStar, public HeGB {

double CalculateCoreMassOnPhase() const { return m_Mass; } // Return m_Mass

double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) { return 0.0; } // Should never be called...
double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate,
const double p_massTransferEfficiencyBeta) { return 0.0; } // Should never be called...

void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch
void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables
Expand Down
6 changes: 4 additions & 2 deletions src/Star.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,8 @@ class Star {
double BindingEnergy_Nanjing() const { return m_Star->BindingEnergy_Nanjing(); }
double BindingEnergy_Kruckow() const { return m_Star->BindingEnergy_Kruckow(); }
double BindingEnergy_Dewi() const { return m_Star->BindingEnergy_Dewi(); }
double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) const { return m_Star->CalculateCriticalMassRatio(p_AccretorIsDegenerate); }
double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate,
const double p_massTransferEfficiencyBeta) const { return m_Star->CalculateCriticalMassRatio(p_AccretorIsDegenerate, p_massTransferEfficiencyBeta); }
double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const { return m_Star->CalculateCriticalMassRatioClaeys14(p_AccretorIsDegenerate); }
double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return m_Star->CalculateCriticalMassRatioHurleyHjellmingWebbink(); }
double CalculateDynamicalTimescale() const { return m_Star->CalculateDynamicalTimescale(); }
Expand Down Expand Up @@ -206,7 +207,8 @@ class Star {

double EvolveOneTimestep(const double p_Dt, const bool p_Force = false);

double InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription) { return m_Star->InterpolateGe20QCrit(p_qCritPrescription); }
double InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription,
const double p_massTransferEfficiencyBeta) { return m_Star->InterpolateGe20QCrit(p_qCritPrescription, p_massTransferEfficiencyBeta); }
void HaltWinds() { m_Star->HaltWinds(); }

void ResolveAccretion(const double p_AccretionMass) { m_Star->ResolveAccretion(p_AccretionMass); }
Expand Down
2 changes: 1 addition & 1 deletion src/constants.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4081,7 +4081,7 @@ const std::map<int, COMPASUnorderedMap<AB_TCoeff, double>> A_COEFF = {
// Note that the radius may contract several times. These points have been removed to facilitate the interpolation, so logR is monotonic.
// Previously, this table contained the commonly used Mass-Radius relation zeta's, however these are functionaly identical to the critical
// mass ratios of the same derivation, and were not used in the code anyway, so they have been replaced by the non-conservative qcrits.
const std::tuple< std::vector<double>, std::vector< std::tuple<std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>>>> GE20_QCRIT_AND_ZETA = {
const std::tuple< std::vector<double>, std::vector< std::tuple<std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>>>> GE20_QCRIT = {
{0.100, 0.130, 0.160, 0.200, 0.220, 0.250, 0.280, 0.320, 0.360, 0.400, 0.450, 0.500, 0.560, 0.630, 0.710, 0.800, 0.890, 1.000, 1.140, 1.300, 1.439, 1.600, 1.800, 2.000, 2.040, 2.500, 3.200, 4.000, 5.000, 6.300, 8.000, 10.000, 13.000, 16.000, 20.000, 25.000, 32.000, 40.000, 50.000, 63.000, 80.000, 100.000},
{
{{-1.1463, -1.1653, -1.1824, -1.1951, -1.1997}, {-0.0027, -0.0018, -0.0021, 0.0020, 0.0033}, {1.5165, 1.5324, 1.5344, 1.5360, 1.5374}, {-94.8767, -109.5290, -103.1992, -390.6250, 9090.9091}, {1.0060, 1.0132, 1.0143, 1.0151, 1.0159}},
Expand Down

0 comments on commit 5e5b8da

Please sign in to comment.