Skip to content

Commit

Permalink
Merge pull request #1291 from TeamCOMPAS/Output
Browse files Browse the repository at this point in the history 
A range of fixes: using core radii for CE survival, fixes to #1255, #1258, #1265
  • Loading branch information
@jeffriley
jeffriley authored Nov 28, 2024
2 parents 7d7b913 + 1f86352 commit e3f130a
Show file tree
Hide file tree
Showing 16 changed files with 205 additions and 147 deletions.
7 changes: 4 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 @@ -847,7 +847,7 @@ Options: { JEANS, ISOTROPIC, CIRCUMBINARY, MACLEOD_LINEAR, ARBITRARY } |br|
Default = ISOTROPIC

**--mass-transfer-fa** |br|
Mass Transfer fraction accreted. |br|
Mass Transfer fraction accreted (beta). |br|
Used when ``--mass-transfer-accretion-efficiency-prescription = FIXED_FRACTION``. |br|
Default = 0.5

Expand Down Expand Up @@ -1122,10 +1122,11 @@ Enable mass loss due to pulsational-pair-instability (PPI). |br|
Default = TRUE

**--pulsational-pair-instability-prescription** |br|
Pulsational pair instability prescription. |br|
Pulsational pair instability prescription (only relevant when using ``--pulsational-pair-instability``). |br|
Options: { HENDRIKS, COMPAS, STARTRACK, MARCHANT, FARMER } |br|
``HENDRIKS`` implements the prescription from Hendriks et al. 2023 |br|
``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, all as implemented in Stevenson et al. 2019. |br|
``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018,
all as implemented in Stevenson et al. 2019. |br|
``FARMER`` follows Farmer et al. 2019 |br|
Default = MARCHANT

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.09.00 Nov 25, 2024**

Improved nuclear timescale mass transfer: the nuclear timescale mass transfer rate is now set by the requirement that the star
ends the time step just filling its Roche lobe.
Fixed several significant mass-transfer issues, such as accretors not gaining mass appropriately and failures
in the root solver for fitting the star into the Roche lobe that were leading to artificial common envelopes and mergers.

**03.08.02 Nov 18, 2024**

Updated implementation of the mass transfer stability critical mass ratio tables from the team of Hongwei Ge.
Expand Down
144 changes: 77 additions & 67 deletions src/BaseBinaryStar.cpp
View file

Large diffs are not rendered by default.

4 changes: 4 additions & 0 deletions src/BaseStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2478,6 +2478,8 @@ double BaseStar::CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu)
*/
double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const {

if (p_Mdot <= 0.0) return 0.0; // nothing to adjust

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

Expand All @@ -2493,6 +2495,8 @@ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(c
logMdotCorrected = logMdotUncorrected;
}

std::cout<<logMdotCorrected<<std::endl;

return PPOW(10.0, logMdotCorrected);
}

Expand Down
5 changes: 5 additions & 0 deletions src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,6 +281,9 @@ class BaseStar {
virtual double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // Default for stars with no convective envelope

virtual double CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const { return 0.0; } // Only defined for MS stars

virtual double CalculateRemnantRadius() const { return Radius(); } // relevant for MS stars, over-written for GB stars


void CalculateSNAnomalies(const double p_Eccentricity);

Expand Down Expand Up @@ -340,6 +343,8 @@ class BaseStar {
void StashSupernovaDetails(const STELLAR_TYPE p_StellarType,
const SSE_SN_RECORD_TYPE p_RecordType = SSE_SN_RECORD_TYPE::DEFAULT) { LOGGING->StashSSESupernovaDetails(this, p_StellarType, p_RecordType); }

virtual double TAMSCoreMass() const { return 0.0; } // except MS stars

virtual void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Default is NO-OP
virtual void UpdateAgeAfterMassLoss() { } // Default is NO-OP

Expand Down
5 changes: 2 additions & 3 deletions src/GiantBranch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,8 @@ class GiantBranch: virtual public BaseStar, public MainSequence {
GiantBranch(){};

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

virtual double CalculateRemnantRadius() const;


double CalculateRemnantRadius() const;

protected:

Expand Down
2 changes: 2 additions & 0 deletions src/HG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#include "HeMS.h"
#include "HeWD.h"


///////////////////////////////////////////////////////////////////////////////////////
// //
// COEFFICIENT AND CONSTANT CALCULATIONS ETC. //
Expand Down Expand Up @@ -861,6 +862,7 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const
}



///////////////////////////////////////////////////////////////////////////////////////
// //
// MASS CALCULATIONS //
Expand Down
80 changes: 39 additions & 41 deletions src/HG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,52 +11,49 @@
#include "GiantBranch.h"


// JR: todo: revisit this one day - sometimes HG works better as GiantBranch, sometimes not...
// Right now it is GiantBranch - figure out which is more appropriate

class BaseStar;
class GiantBranch;

class HG: virtual public BaseStar, public GiantBranch {

public:

HG() { m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; };

HG(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), GiantBranch(p_BaseStar) {
m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; // Set stellar type
m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; // Set stellar type
if (p_Initialise) Initialise(); // Initialise if required
}

HG* Clone(const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
HG* clone = new HG(*this, p_Initialise);
clone->SetPersistence(p_Persistence);
return clone;
HG* clone = new HG(*this, p_Initialise);
clone->SetPersistence(p_Persistence);
return clone;
}

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;
HG* clone = new HG(p_Star, p_Initialise);
clone->SetPersistence(p_Persistence);
return clone;
}

MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::B; } // Always case B


protected:

void Initialise() {

m_Tau = 0.0; // Start of phase

// update stellar properties at start of HG phase (since core definition changes)
CalculateGBParams();
CalculateTimescales();
// Initialise timescales
// Initialise timescales
m_Age = m_Timescales[static_cast<int>(TIMESCALE::tMS)]; // Set age appropriately

// update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass
// (only relevant if RetainCoreMassDuringCaseAMassTransfer())
// (only relevant if RetainCoreMassDuringCaseAMassTransfer())
if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) {
double desiredCoreMass = std::min(m_Mass, MinimumCoreMass()); // desired core mass
m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass); // use root finder to find new core mass estimate
Expand All @@ -69,56 +66,56 @@ class HG: virtual public BaseStar, public GiantBranch {
}
EvolveOnPhase(0.0);
}


// member functions - alphabetically
double CalculateCOCoreMassAtPhaseEnd() const { return 0.0; } // McCO(HG) = 0.0
double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(HG) = 0.0

double CalculateCoreMassAt2ndDredgeUp(const DBL_VECTOR &p_GBParams) { return p_GBParams[static_cast<int>(GBP::McDU)]; } // NO-OP
double CalculateCoreMassAtPhaseEnd(const double p_Mass) const;
double CalculateCoreMassAtPhaseEnd() const { return CalculateCoreMassAtPhaseEnd(m_Mass0); } // Use class member variables
double CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) const;
double CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // Use class member variables
double CalculateCoreMassOnPhaseIgnoringPreviousCoreMass(const double p_Mass, const double p_Time) const; // Ignore previous core mass constraint when computing expected core mass

double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const;
double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.25; } // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.

double CalculateHeCoreMassAtPhaseEnd() const { return m_CoreMass; } // McHe(HG) = Core Mass
double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // McHe(HG) = Core Mass

double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 1.0 - m_Metallicity; }
double CalculateHeliumAbundanceCoreOnPhase() const { return 1.0 - m_Metallicity; } // Use class member variables
double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 1.0 - m_Metallicity; }
double CalculateHeliumAbundanceCoreOnPhase() const { return 1.0 - m_Metallicity; } // Use class member variables

double CalculateHeliumAbundanceSurfaceAtPhaseEnd() const { return CalculateHeliumAbundanceSurfaceOnPhase(); }
double CalculateHeliumAbundanceSurfaceOnPhase() const { return m_InitialHeliumAbundance; } // Use class member variables
double CalculateHeliumAbundanceSurfaceOnPhase() const { return m_InitialHeliumAbundance; } // Use class member variables

double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return CalculateHydrogenAbundanceCoreOnPhase(); }
double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const;
double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; } // Star has exhausted hydrogen in its core
double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return CalculateHydrogenAbundanceCoreOnPhase(); }
double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const;
double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; } // Star has exhausted hydrogen in its core

double CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const { return CalculateHydrogenAbundanceSurfaceOnPhase(); }
double CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const { return CalculateHydrogenAbundanceSurfaceOnPhase(); }
double CalculateHydrogenAbundanceSurfaceOnPhase() const { return m_InitialHydrogenAbundance; } // Use class member variables



double CalculateLambdaDewi() const;
double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const;
double CalculateLambdaNanjingEnhanced(const int p_MassIndex, const STELLAR_POPULATION p_StellarPop) const;

double CalculateLuminosityAtPhaseEnd(const double p_Mass) const;
double CalculateLuminosityAtPhaseEnd() const { return CalculateLuminosityAtPhaseEnd(m_Mass0);} // Use class member variables
double CalculateLuminosityOnPhase(const double p_Age, const double p_Mass) const;
double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age, m_Mass0); } // Use class member variables

double CalculateMassTransferRejuvenationFactor() { return 1.0; }

double CalculateRadiusAtPhaseEnd(const double p_Mass) const;
double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass); } // Use class member variables
double CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const double p_RZAMS) const;
double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass0, m_Tau, m_RZAMS0); } // Use class member variables

double CalculateRho(const double p_Mass) const;

double CalculateTauAtPhaseEnd() const { return 1.0; } // tau = 1.0 at end of HG
Expand All @@ -141,6 +138,7 @@ class HG: virtual public BaseStar, public GiantBranch {
bool ShouldEvolveOnPhase() const { return (utils::Compare(m_Age, m_Timescales[static_cast<int>(TIMESCALE::tBGB)]) < 0); } // Evolve on HG phase if age < Base Giant Branch timescale
bool ShouldSkipPhase() const { return false; } // Never skip HG phase

double TAMSCoreMass() const { return 0.0; }
void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Nothing to do for stars beyond the Main Sequence for now
void UpdateAgeAfterMassLoss(); // Per Hurley et al. 2000, section 7.1

Expand Down
16 changes: 16 additions & 0 deletions src/HeHG.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "HeHG.h"
#include "HeGB.h"
#include "COWD.h"
#include "HeWD.h"


/*
Expand Down Expand Up @@ -234,6 +235,21 @@ double HeHG::CalculatePerturbationMu() const {
return std::max(5.0 * ((McMax - m_CoreMass) / McMax), 0.0); //return non-negative value to avoid round-off issues
}

/*
* Calculate radius of the remnant the star would become if it lost all of its
* envelope immediately (i.e. M = Mc, coreMass)
*
* Hurley et al. 2000, end of section 6
*
*
* double CalculateRemnantRadius()
*
* @return Radius of remnant core in Rsol
*/
double HeHG::CalculateRemnantRadius() const {
return HeWD::CalculateRadiusOnPhase_Static(m_CoreMass);
}


///////////////////////////////////////////////////////////////////////////////////////
// //
Expand Down
1 change: 1 addition & 0 deletions src/HeHG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ class HeHG: virtual public BaseStar, public HeMS {
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);

double CalculateRemnantRadius() const;

protected:

Expand Down
2 changes: 2 additions & 0 deletions src/HeMS.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,8 @@ class HeMS: virtual public BaseStar, public TPAGB {
static DBL_DBL CalculateRadiusAtPhaseEnd_Static(const double p_Mass, const double p_Luminosity);
static double CalculateRadiusAtZAMS_Static(const double p_Mass);
static double CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Tau);

double CalculateRemnantRadius() const { return Radius(); }

MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::OTHER; } // Not A, B, C, or NONE

Expand Down
Loading

0 comments on commit e3f130a

Please sign in to comment.