Skip to content

Commit

Permalink
Merge pull request #1150 from TeamCOMPAS/MassTransferODE
Browse files Browse the repository at this point in the history 
Mass transfer ODE
  • Loading branch information
@jeffriley
jeffriley authored Jun 22, 2024
2 parents 723e278 + 944205d commit 1766d94
Show file tree
Hide file tree
Showing 5 changed files with 46 additions and 32 deletions.
54 changes: 30 additions & 24 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1810,19 +1810,18 @@ double BaseBinaryStar::CalculateRocheLobeRadius_Static(const double p_MassPrimar
* which is the fraction of specific angular momentum with which the non-accreted mass leaves the system.
* Macleod_linear comes from Willcox et al. (2022)
*
* Updates class member variable m_Error JR: todo: revisit error handling (this could be a const function)
*
*
* Calculation is based on user-specified Angular Momentum Loss prescription
*
*
* double CalculateGammaAngularMomentumLoss(const double p_DonorMass, const double p_AccretorMass)
* double CalculateGammaAngularMomentumLoss_Static(const double p_DonorMass, const double p_AccretorMass, const bool p_IsAccretorDegenerate)
*
* @param [IN] p_DonorMass The mass of the donor (Msol)
* @param [IN] p_AccretorMass The mass of the accretor (Msol)
* @param [IN] p_IsAccretorDegenerate True if the accretor is a degenerate star, false otherwise (need to know up front to keep this function static)
* @return The fraction of specific angular momentum with which the non-accreted mass leaves the system
*/
double BaseBinaryStar::CalculateGammaAngularMomentumLoss(const double p_DonorMass, const double p_AccretorMass) {
double BaseBinaryStar::CalculateGammaAngularMomentumLoss_Static(const double p_DonorMass, const double p_AccretorMass, const bool p_IsAccretorDegenerate) {

double gamma;

Expand All @@ -1835,8 +1834,8 @@ double BaseBinaryStar::CalculateGammaAngularMomentumLoss(const double p_DonorMas
// interpolate in separation between a_acc and a_L2, both normalized to units of separation a
double aL2 = std::sqrt(M_SQRT2); // roughly, coincides with CIRCUMBINARY_RING def above
double aAcc = 1.0 / (1.0 + q);
double fMacleod = m_Accretor->IsDegenerate()
? OPTIONS->MassTransferJlossMacLeodLinearFractionDegen()
double fMacleod = p_IsAccretorDegenerate
? OPTIONS->MassTransferJlossMacLeodLinearFractionDegen()
: OPTIONS->MassTransferJlossMacLeodLinearFractionNonDegen();
double aGamma = aAcc + (aL2 - aAcc)*fMacleod;
gamma = aGamma * aGamma * (1.0 + q) * (1.0 + q) / q;
Expand All @@ -1845,8 +1844,6 @@ double BaseBinaryStar::CalculateGammaAngularMomentumLoss(const double p_DonorMas
case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::ARBITRARY : gamma = OPTIONS->MassTransferJloss(); break;
default: // unknown mass transfer angular momentum loss prescription - shouldn't happen
gamma = 1.0; // default value
m_Error = ERROR::UNKNOWN_MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION; // set error value
SHOW_WARN(m_Error); // warn that an error occurred
}

return gamma;
Expand Down Expand Up @@ -1875,25 +1872,34 @@ double BaseBinaryStar::CalculateMassTransferOrbit(const double p
BinaryConstituentStar& p_Accretor,
const double p_FractionAccreted) {

double semiMajorAxis = m_SemiMajorAxis; // new semi-major axis value - default is no change
double massA = p_Accretor.Mass(); // accretor mass
double massD = p_DonorMass; // donor mass
double jLoss; // specific angular momentum carried away by non-conservative mass transfer

double semiMajorAxis = m_SemiMajorAxis;

if (utils::Compare(p_DeltaMassDonor, 0.0) < 0) { // mass loss from donor?
// yes
int numberIterations = fmax( floor (fabs(p_DeltaMassDonor/(MAXIMUM_MASS_TRANSFER_FRACTION_PER_STEP * massD))), 1.0); // number of iterations
double dM = p_DeltaMassDonor / numberIterations; // mass change per time step

for(int i = 0; i < numberIterations ; i++) {
jLoss = CalculateGammaAngularMomentumLoss(massD, massA);
semiMajorAxis = semiMajorAxis + (((-2.0 * dM / massD) * (1.0 - (p_FractionAccreted * (massD / massA)) - ((1.0 - p_FractionAccreted) * (jLoss + 0.5) * (massD / (massA + massD))))) * semiMajorAxis);

massD = massD + dM;
massA = massA - (dM * p_FractionAccreted);
}
}
controlled_stepper_type controlled_stepper;
state_type x(1);
x[0] = semiMajorAxis;

// Use boost adaptive ODE solver for speed and accuracy
struct ode
{
double p_MassDonor0, p_MassAccretor0, p_FractionAccreted;
bool p_IsAccretorDegenerate;
ode(double massDonor0, double massAccretor0, double fractionAccreted, bool isAccretorDegenerate) : p_MassDonor0(massDonor0), p_MassAccretor0(massAccretor0), p_FractionAccreted(fractionAccreted), p_IsAccretorDegenerate(isAccretorDegenerate) {}

void operator()( state_type const& x , state_type& dxdt , double p_MassChange ) const
{
double massD = p_MassDonor0 + p_MassChange;
double massA = p_MassAccretor0 - p_MassChange * p_FractionAccreted;
double jLoss = CalculateGammaAngularMomentumLoss_Static(massD, massA, p_IsAccretorDegenerate);
dxdt[0] = (-2.0 / massD) * (1.0 - (p_FractionAccreted * (massD / massA)) - ((1.0 - p_FractionAccreted) * (jLoss + 0.5) * (massD / (massA + massD)))) * x[0];
}
};

integrate_adaptive( controlled_stepper , ode { p_DonorMass, p_Accretor.Mass(), p_FractionAccreted, m_Accretor->IsDegenerate() }, x , 0.0 , p_DeltaMassDonor , p_DeltaMassDonor/1000.0 );
semiMajorAxis = x[0];
}

return semiMajorAxis;
}

Expand Down
10 changes: 6 additions & 4 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "BinaryConstituentStar.h"

#include <boost/math/tools/roots.hpp>
#include <boost/numeric/odeint.hpp>

//#include <boost/math/special_functions/next.hpp> // For float_distance.
//#include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.
Expand Down Expand Up @@ -425,8 +426,9 @@ class BaseBinaryStar {
double CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity);
double CalculateDOmegaTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, const double p_M1, const double p_R1, const double p_I1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity);
double CalculateDSemiMajorAxisTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity);

double CalculateGammaAngularMomentumLoss(const double p_DonorMass, const double p_AccretorMass);

static double CalculateGammaAngularMomentumLoss_Static(const double p_DonorMass, const double p_AccretorMass, const bool p_IsAccretorDegenerate);
double CalculateGammaAngularMomentumLoss(const double p_DonorMass, const double p_AccretorMass) { return CalculateGammaAngularMomentumLoss_Static(p_DonorMass, p_AccretorMass, m_Accretor->IsDegenerate()); }
double CalculateGammaAngularMomentumLoss() { return CalculateGammaAngularMomentumLoss(m_Donor->Mass(), m_Accretor->Mass()); }


Expand All @@ -449,7 +451,7 @@ class BaseBinaryStar {
const double p_Mass,
const double p_SemiMajorAxis) const { return -(G_AU_Msol_yr * p_Mu * p_Mass) / (2.0 * p_SemiMajorAxis); }

double CalculateZetaRocheLobe(const double p_jLoss, const double p_beta) const;
double CalculateZetaRocheLobe(const double p_jLoss, const double p_beta) const;

double CalculateTimeToCoalescence(double a0, double e0, double m1, double m2) const;

Expand Down Expand Up @@ -768,7 +770,7 @@ class BaseBinaryStar {

return root.first + (root.second - root.first) / 2.0; // midway between brackets (could return brackets...)
}

};

#endif // __BaseBinaryStar_h__
4 changes: 3 additions & 1 deletion src/changelog.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1216,9 +1216,11 @@
// - Reduced MAXIMUM_MASS_TRANSFER_FRACTION_PER_STEP to 0.0001 to improve accuracy of orbital separation updates following mass transfer
// - Corrected temperature units in Picker formula for Tonset used in the calculation of the convective envelope mass
// - Code cleanup
// 02.49.05 IM - June 22, 2024 - Enhancement:
// - Replaced fixed-step, first-order integrator for orbital change after mass transfer with an adaptive-step, higher-order ODE integrator for improved speed and accuracy


const std::string VERSION_STRING = "02.49.04";
const std::string VERSION_STRING = "02.49.05";



Expand Down
3 changes: 0 additions & 3 deletions src/constants.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -326,9 +326,6 @@ constexpr int MAX_TIMESTEP_RETRIES = 30;
constexpr double MAXIMUM_MASS_LOSS_FRACTION = 0.01; // Maximum allowable mass loss - 1.0% (of mass) expressed as a fraction
constexpr double MAXIMUM_RADIAL_CHANGE = 0.01; // Maximum allowable radial change - 1% (of radius) expressed as a fraction
constexpr double MINIMUM_MASS_SECONDARY = 4.0; // Minimum mass of secondary to evolve

constexpr double MAXIMUM_MASS_TRANSFER_FRACTION_PER_STEP= 0.0001; // Maximal fraction of donor mass that can be transferred in one step of stable mass transfer

constexpr double LAMBDA_NANJING_ZLIMIT = 0.0105; // Metallicity cutoff for Nanjing lambda calculations
constexpr double LAMBDA_NANJING_POPI_Z = 0.02; // Population I metallicity in Xu & Li (2010)
constexpr double LAMBDA_NANJING_POPII_Z = 0.001; // Population II metallicity in Xu & Li (2010)
Expand Down
7 changes: 7 additions & 0 deletions src/typedefs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
#define __typedefs_h__

#include "constants.h"
#include <boost/math/tools/roots.hpp>
#include <boost/numeric/odeint.hpp>

// JR: todo: clean this up and document it better

Expand Down Expand Up @@ -306,4 +308,9 @@ typedef struct StellarCEDetails { // Common Envelope detail
} StellarCEDetailsT; // was CommonEnvelopeDetailsT;


// For boost ODE integrators, see https://www.boost.org/doc/libs/1_83_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial/harmonic_oscillator.html
typedef std::vector< double > state_type;
typedef boost::numeric::odeint::runge_kutta_cash_karp54< state_type > error_stepper_type;
typedef boost::numeric::odeint::controlled_runge_kutta< error_stepper_type > controlled_stepper_type;

#endif // __typedefs_h__

0 comments on commit 1766d94

Please sign in to comment.