Skip to content

Commit

Permalink
Merge pull request #1230 from TeamCOMPAS/PickerFixes
Browse files Browse the repository at this point in the history 
Multiple changes to the prescriptions for the masses and radii of convective and radiative layers in stars
  • Loading branch information
@ilyamandel
ilyamandel authored Sep 30, 2024
2 parents aaf157f + 9844fa9 commit 66cfa99
Show file tree
Hide file tree
Showing 11 changed files with 143 additions and 51 deletions.
2 changes: 1 addition & 1 deletion src/EAGB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -855,7 +855,7 @@ double EAGB::CalculateRemnantRadius() const {
double R1, R2;
std::tie(R1, R2) = HeGB::CalculateRadiusOnPhase_Static(m_HeCoreMass, CalculateRemnantLuminosity());

return std::min(R1, R2);
return std::min(m_Radius, std::min(R1, R2));
}


Expand Down
58 changes: 41 additions & 17 deletions src/GiantBranch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -679,6 +679,32 @@ double GiantBranch::CalculateRemnantRadius() const {
#undef massCutoffs
}

/*
* Calculate the radial extent of the convective outer envelope
*
* Combination of Hurley et al. 2000, end of sec. 7.2, and Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 39-40
*
*
* double CalculateRadialExtentConvectiveEnvelope()
*
* @return Radial extent of the convective outer envelope
*/
double GiantBranch::CalculateRadialExtentConvectiveEnvelope() const{
double convectiveEnvelopeMass, convectiveEnvelopeMassMax;
std::tie(convectiveEnvelopeMass, convectiveEnvelopeMassMax) = CalculateConvectiveEnvelopeMass();
if (utils::Compare(convectiveEnvelopeMass, 0.0) <= 0 || utils::Compare(convectiveEnvelopeMassMax, 0.0) <= 0 ) return 0.0; // massless convective envelope has zero radial extent

double convectiveCoreMass = CalculateConvectiveCoreMass();
double convectiveCoreRadius = CalculateConvectiveCoreRadius();
// assume that the final radiative intershell (if any) would have a density that is a geometric mean of the core density and total density
double radiativeIntershellMass = m_Mass - convectiveCoreMass - convectiveEnvelopeMassMax;
double convectiveCoreRadiusCubed = convectiveCoreRadius * convectiveCoreRadius * convectiveCoreRadius;
double radiativeIntershellDensity = 1.0 / (4.0 /3.0 * M_PI) * std::sqrt(convectiveCoreMass / convectiveCoreRadiusCubed * m_Mass / m_Radius / m_Radius / m_Radius);
double outerEdgeRadiativeIntershellCubed = radiativeIntershellMass / (4.0 / 3.0 * M_PI * radiativeIntershellDensity) + convectiveCoreRadiusCubed;

return std::sqrt(convectiveEnvelopeMass/convectiveEnvelopeMassMax) * (m_Radius - std::cbrt(outerEdgeRadiativeIntershellCubed));
}


///////////////////////////////////////////////////////////////////////////////////////
// //
Expand All @@ -700,9 +726,7 @@ double GiantBranch::CalculateRemnantRadius() const {
*/
double GiantBranch::CalculateCoreMassAtBAGB(const double p_Mass) const {
#define b m_BnCoefficients // for convenience and readability - undefined at end of function

return std::sqrt(std::sqrt((b[36] * PPOW(p_Mass, b[37])) + b[38])); // sqrt() is much faster than PPOW()

#undef b
}

Expand All @@ -723,9 +747,7 @@ double GiantBranch::CalculateCoreMassAtBAGB(const double p_Mass) const {
*/
double GiantBranch::CalculateCoreMassAtBAGB_Static(const double p_Mass, const DBL_VECTOR &p_BnCoefficients) {
#define b p_BnCoefficients // for convenience and readability - undefined at end of function

return std::sqrt(std::sqrt((b[36] * PPOW(p_Mass, b[37])) + b[38])); // sqrt() is much faster than PPOW()

#undef b
}

Expand Down Expand Up @@ -1049,9 +1071,7 @@ double GiantBranch::CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPre
*/
DBL_DBL GiantBranch::CalculateConvectiveEnvelopeMass() const {

double MinterfMcoref = -0.023 * m_Log10Metallicity - 0.0023; // Eq. (8) of Picker+ 2024
double Tonset = -129.7 * m_Log10Metallicity * m_Log10Metallicity - 920.1 * m_Log10Metallicity + 2887.1; // Eq. (6) of Picker+ 2024
Tonset /= TSOL; // convert to solar units
double MinterfMcoref = -0.023 * m_Log10Metallicity - 0.0023; // eq. (8) of Picker+ 2024

// We need the temperature of the star just after BAGB, which is the temperature at the
// start of the EAGB phase. Since we are on the giant branch here, we can clone this
Expand All @@ -1065,21 +1085,25 @@ DBL_DBL GiantBranch::CalculateConvectiveEnvelopeMass() const {
// underlying object is, we cast it to EAGB&.

EAGB *clone = EAGB::Clone(static_cast<EAGB&>(const_cast<GiantBranch&>(*this)), OBJECT_PERSISTENCE::EPHEMERAL);
clone->UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true); // Otherwise, temperature not updated
clone->UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true); // otherwise, temperature not updated
double Tmin = clone->Temperature(); // get temperature of clone
delete clone; clone = nullptr; // return the memory allocated for the clone

// Assume Tonset is always 0.1 dex hotter than Tmin, to avoid issues caused by differences between temperatures
// in MESA models (used in Picker+ fits) and Pols models (used in Hurley+ SSE tracks), rather than Eq. (6) of Picker+ 2024
double Tonset = 1.2589 * Tmin;

double mCoreFinal = CalculateCoreMassAtBAGB(m_Mass0);
double mConvMax = std::max(m_Mass - mCoreFinal * (1.0 + MinterfMcoref), 0.0); // eq. (9) of Picker+ 2024

double McoreFinal = CalculateCoreMassAtBAGB(m_Mass0);
double MconvMax = std::max(m_Mass - McoreFinal * (1.0 + MinterfMcoref), 0.0); // Eq. (9) of Picker+ 2024
if( utils::Compare(McoreFinal,1.5) < 0 ) // Picker+ 2024 fits were only made for stars above 8.0 solar masses, with runs down to 5.0 solar masses, so using the final core mass as an approximate threshold of validity
MconvMax = std::max(m_Mass - McoreFinal, 0.0); // unlike massive stars, intermediate-mass stars have almost no radiative intershell at maximum convective envelope extent
double convectiveEnvelopeMass = MconvMax / (1.0 + exp(4.6 * (Tmin + Tonset - 2.0 * m_Temperature) / (Tmin - Tonset))); // Eq. (7) of Picker+ 2024
// Picker+ 2024 fits were only made for stars above 8.0 solar masses, with runs down to 5.0 solar masses,
// so using the final core mass as an approximate threshold of validity
if(utils::Compare(mCoreFinal, 1.5) < 0) mConvMax = std::max(m_Mass - mCoreFinal, 0.0); // unlike massive stars, intermediate-mass stars have almost no radiative intershell at maximum convective envelope extent

if(utils::Compare(Temperature(), Tmin) <= 0)
convectiveEnvelopeMass = MconvMax; // If already on the AGB, the convective envelope mass should be set to its maximum value
double convectiveEnvelopeMass = mConvMax / (1.0 + exp(4.6 * (Tmin + Tonset - 2.0 * m_Temperature) / (Tmin - Tonset))); // eq. (7) of Picker+ 2024

return std::tuple<double, double> (convectiveEnvelopeMass, MconvMax);
} // /*ILYA*/ check consistency with HG convective envelope radii and masses from Hurley+ 2002, 2000
return std::tuple<double, double> (convectiveEnvelopeMass, mConvMax);
}


///////////////////////////////////////////////////////////////////////////////////////
Expand Down
16 changes: 6 additions & 10 deletions src/GiantBranch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,15 @@ class GiantBranch: virtual public BaseStar, public MainSequence {

GiantBranch(const BaseStar &p_BaseStar) : BaseStar(p_BaseStar), MainSequence(p_BaseStar) {}

virtual double CalculateRemnantRadius() const;


protected:


// member functions - alphabetically (sort of - some are grouped by functionality)
double CalculateConvectiveCoreMass() const { return m_CoreMass; }
double CalculateConvectiveCoreRadius () const { return std::min(CalculateRemnantRadius (), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000
double CalculateConvectiveCoreMass() const { return m_CoreMass; }
double CalculateConvectiveCoreRadius () const { return std::min(CalculateRemnantRadius (), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000
DBL_DBL CalculateConvectiveEnvelopeMass() const;
static double CalculateCoreMassAt2ndDredgeUp_Static(const double p_McBAGB);
double CalculateCoreMassAtBAGB(const double p_Mass) const;
Expand All @@ -49,7 +51,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence {

void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams);
static void CalculateGBParams_Static(const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_BnCoefficients, DBL_VECTOR &p_GBParams);
void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables
void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables

static double CalculateHRateConstant_Static(const double p_Mass);

Expand Down Expand Up @@ -92,12 +94,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence {

double CalculatePerturbationMu() const;

double CalculateRadialExtentConvectiveEnvelope() const {
// combination of Hurley et al. 2000, end of sec. 7.2, and Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 39-40
double envMass, envMassMax;
std::tie(envMass, envMassMax) = CalculateConvectiveEnvelopeMass();
return (std::sqrt(envMass / envMassMax) * (m_Radius - CalculateConvectiveCoreRadius()));
}
double CalculateRadialExtentConvectiveEnvelope() const;

double CalculateRadiusAtHeIgnition(const double p_Mass) const;
double CalculateRadiusOnPhase(const double p_Mass, const double p_Luminosity) const { return CalculateRadiusOnPhase_Static(p_Mass, p_Luminosity, m_BnCoefficients); }
Expand All @@ -114,7 +111,6 @@ class GiantBranch: virtual public BaseStar, public MainSequence {
virtual double CalculateRemnantLuminosity() const;
STELLAR_TYPE CalculateRemnantTypeByMuller2016(const double p_COCoreMass);

virtual double CalculateRemnantRadius() const;

double CalculateThermalMassLossRate() const { return (m_Mass - m_CoreMass) / CalculateThermalTimescale(); } // Use class member variables

Expand Down
2 changes: 1 addition & 1 deletion src/HG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ class HG: virtual public BaseStar, public GiantBranch {
return clone;
}

static HG* Clone(HG p_Star, const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
static HG* Clone(HG& p_Star, const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
HG* clone = new HG(p_Star, p_Initialise);
clone->SetPersistence(p_Persistence);
return clone;
Expand Down
1 change: 0 additions & 1 deletion src/HeGB.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,6 @@ class HeGB: virtual public BaseStar, public HeHG {

static double CalculateLuminosityOnPhase_Static(const double p_CoreMass, const double p_GBPB, const double p_GBPD);

double CalculateRadialExtentConvectiveEnvelope() const { return FGB::CalculateRadialExtentConvectiveEnvelope(); }
static DBL_DBL CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Luminosity);


Expand Down
9 changes: 3 additions & 6 deletions src/HeHG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,14 +37,15 @@ class HeHG: virtual public BaseStar, public HeMS {


// member functions
double CalculateConvectiveCoreRadius () const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000
static void CalculateGBParams_Static(const double p_Mass0, const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_BnCoefficients, DBL_VECTOR &p_GBParams);


protected:

void Initialise() {
m_Tau = 0.0; // Start of phase
CalculateTimescales(); // Initialise timescales
CalculateTimescales(); // Initialise timescales
// Age for HeHG is calculated before switching -
// can get here via EvolveOneTimestep() and ResolveEnvelopeLoss(),
// and Age is calculated differently in those cases
Expand All @@ -67,8 +68,6 @@ class HeHG: virtual public BaseStar, public HeMS {

double CalculateConvectiveCoreMass() const { return m_CoreMass; }

double CalculateConvectiveCoreRadius () const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000

double CalculateCoreMassAtBAGB() const { return m_Mass0; } // McBAGB = M0 (Hurely et al. 2000, discussion just before eq 89)
double CalculateCoreMassAtPhaseEnd() const { return m_CoreMass; } // NO-OP
double CalculateCoreMassOnPhase() const { return m_COCoreMass; } // Mc(HeMS) = McCOMass
Expand All @@ -95,8 +94,6 @@ class HeHG: virtual public BaseStar, public HeMS {
double CalculatePerturbationMu() const;
double CalculatePerturbationMuAtPhaseEnd() const { return m_Mu; } // NO-OP

double CalculateRadialExtentConvectiveEnvelope() const { return HG::CalculateRadialExtentConvectiveEnvelope(); }

double CalculateRadiusAtPhaseEnd() const { return m_Radius; } // NO-OP
double CalculateRadiusOnPhase() const;

Expand Down
57 changes: 57 additions & 0 deletions src/HeMS.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,34 @@ double HeMS::CalculateLuminosityOnPhase_Static(const double p_Mass, const double
// //
///////////////////////////////////////////////////////////////////////////////////////

/*
* Calculate convective core radius
*
* Assume equal to total radius at start (for continuity with stripped CHeB or HG star), continuity with HeHG at end of phase, linear growth
*
*
* double CalculateConvectiveCoreRadius()
*
* @return Convective core radius (solar radii)
*/
double HeMS::CalculateConvectiveCoreRadius() const {

// We need core radius at end of phase, which is just the core radius at the start of the HeHG phase.
// Since we are on the He main sequence here, we can clone this object as an HeHG object
// and, as long as it is initialised (to correctly set Tau to 0.0 on the HeHG phase),
// we can query the cloned object for its core mass.
//
// The clone should not evolve, and so should not log anything, but to be sure the
// clone does not participate in logging, we set its persistence to EPHEMERAL.

HeHG *clone = HeHG::Clone(static_cast<HeHG&>(const_cast<HeMS&>(*this)), OBJECT_PERSISTENCE::EPHEMERAL);
double finalConvectiveCoreRadius = clone->CalculateConvectiveCoreRadius(); // get core radius from clone
delete clone; clone = nullptr; // return the memory allocated for the clone

double initialConvectiveCoreRadius = m_Radius;
return (initialConvectiveCoreRadius - m_Tau * (initialConvectiveCoreRadius - finalConvectiveCoreRadius));
}


/*
* Calculate radius at ZAMS for a Helium Main Sequence star
Expand Down Expand Up @@ -227,6 +255,35 @@ double HeMS::CalculateRadiusAtPhaseEnd_Static(const double p_Mass) {
///////////////////////////////////////////////////////////////////////////////////////


/*
* Calculate convective core mass
*
* Assume equal to total mass at start (for continuity with stripped CHeB or HG star), continuity with HeHG at end of phase, linear growth
*
*
* double CalculateConvectiveCoreMass()
*
* @return Convective core mass (solar masses)
*/
double HeMS::CalculateConvectiveCoreMass() const {

// We need core mass at end of phase, which is just the core mass at the start of the HeHG phase.
// Since we are on the He main sequence here, we can clone this object as an HeHG object
// and, as long as it is initialised (to correctly set Tau to 0.0 on the HeHG phase),
// we can query the cloned object for its core mass.
//
// The clone should not evolve, and so should not log anything, but to be sure the
// clone does not participate in logging, we set its persistence to EPHEMERAL.

HeHG *clone = HeHG::Clone(static_cast<HeHG&>(const_cast<HeMS&>(*this)), OBJECT_PERSISTENCE::EPHEMERAL, true);
double finalConvectiveCoreMass = clone->CoreMass(); // get core mass from clone
delete clone; clone = nullptr; // return the memory allocated for the clone

double initialConvectiveCoreMass = m_Mass;
return (initialConvectiveCoreMass - m_Tau * (initialConvectiveCoreMass - finalConvectiveCoreMass));
}


/*
* Calculate rejuvenation factor for stellar age based on mass lost/gained during mass transfer
*
Expand Down
Loading

0 comments on commit 66cfa99

Please sign in to comment.