Merge pull request #1282 from TeamCOMPAS/BugFix

Multiple rotation-related fixes to 03.08.01 (units, initialisation, min->max typo, retain omega on envelope loss)

online-docs/pages/whats-new.rst

@@ -3,6 +3,19 @@ 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.00 Nov 18, 2024**
+
+Improved the treatment of stellar rotation (with further corrections in 03.08.01):
+
+* Assume rigid body rotation
+* Keep the angular moment of a star constant when there is no mass loss
+* When a star with radius r and angular frequency omega loses mass dm through winds or mass transfer, it loses angular momentum dL = 
+(2/3) dm r^2 omega
+* (However, angular momentum never drops below zero)
+* When a star loses its envelope, the remaining core is assumed to rotate with the same rotation rate as the preceding star
+* When a star of mass m and radius r gains mass dm through accretion, it gain angular momentum dL = dm \sqrt{G m r}
+* If initial binary rotation is fast enough for a star to be CHE, it is set to that rotation frequency without regard for the tidal 
+prescription; CHE stars remain tidally locked if the tidal prescription is NONE
 
 **03.07.01 Oct 23, 2024**
 

src/BaseBinaryStar.cpp

@@ -2417,7 +2417,7 @@ void BaseBinaryStar::ResolveMassChanges() {
             // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius
             double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ?
                 massChange * sqrt(G_AU_Msol_yr * m_Star1->Mass() * m_Star1->Radius() * RSOL_TO_AU) :
-                (2.0 / 3.0) * massChange * m_Star1->Radius() * m_Star1->Radius() * m_Star1->Omega();
+                (2.0 / 3.0) * massChange * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Omega();
             // update mass of star according to mass loss and mass transfer, then update age accordingly
             (void)m_Star1->UpdateAttributes(massChange, 0.0);                                           // update mass for star
             m_Star1->UpdateInitialMass();                                                               // update effective initial mass of star (MS, HG & HeMS)
@@ -2440,7 +2440,7 @@ void BaseBinaryStar::ResolveMassChanges() {
             // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius
             double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ?
                 massChange * sqrt(G_AU_Msol_yr * m_Star2->Mass() * m_Star2->Radius() * RSOL_TO_AU) :
-                (2.0 / 3.0) * massChange * m_Star2->Radius() * m_Star2->Radius() * m_Star2->Omega();
+                (2.0 / 3.0) * massChange * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Omega();
             // update mass of star according to mass loss and mass transfer, then update age accordingly
             (void)m_Star2->UpdateAttributes(massChange, 0.0);                                           // update mass for star
             m_Star2->UpdateInitialMass();                                                               // update effective initial mass of star (MS, HG & HeMS)

src/BaseStar.cpp

@@ -111,13 +111,6 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
     m_LZAMS                                    = CalculateLuminosityAtZAMS(m_MZAMS);
     m_TZAMS                                    = CalculateTemperatureOnPhase_Static(m_LZAMS, m_RZAMS);
 
-    m_OmegaCHE                                 = CalculateOmegaCHE(m_MZAMS, m_Metallicity);
-
-    m_OmegaZAMS                                = p_RotationalFrequency >= 0.0                           // valid rotational frequency passed in?
-                                                    ? p_RotationalFrequency                             // yes - use it
-                                                    : CalculateZAMSAngularFrequency(m_MZAMS, m_RZAMS);  // no - calculate it
-    m_AngularMomentum                          = CalculateMomentOfInertiaAU() * m_OmegaZAMS;
-
     // Initial abundances
     m_InitialHeliumAbundance                   = CalculateInitialHeliumAbundance();
     m_HeliumAbundanceCore                      = m_InitialHeliumAbundance;
@@ -143,6 +136,12 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
     m_Radius                                   = m_RZAMS;
     m_Temperature                              = m_TZAMS;
 	m_ComponentVelocity						   = Vector3d();
+
+    m_OmegaCHE                                 = CalculateOmegaCHE(m_MZAMS, m_Metallicity);
+    m_OmegaZAMS                                = p_RotationalFrequency >= 0.0                           // valid rotational frequency passed in?
+                                                    ? p_RotationalFrequency                             // yes - use it
+                                                    : CalculateZAMSAngularFrequency(m_MZAMS, m_RZAMS);  // no - calculate it
+    m_AngularMomentum                          = CalculateMomentOfInertiaAU() * m_OmegaZAMS;
 
     m_CoreMass                                 = DEFAULT_INITIAL_DOUBLE_VALUE;
     m_COCoreMass                               = DEFAULT_INITIAL_DOUBLE_VALUE;
@@ -2970,8 +2969,8 @@ void BaseStar::ResolveMassLoss(const bool p_UpdateMDt) {
 
         double mass = CalculateMassLossValues(true, p_UpdateMDt);                                   // calculate new values assuming mass loss applied
 
-        double angularMomentumChange = (2.0/3.0) * (mass - m_Mass) * m_Radius * m_Radius * Omega();
-
+        double angularMomentumChange = (2.0/3.0) * (mass - m_Mass) * m_Radius * RSOL_TO_AU * m_Radius * RSOL_TO_AU * Omega();
+                
         // JR: this is here to keep attributes in sync BSE vs SSE
         // Supernovae are caught in UpdateAttributesAndAgeOneTimestep() (hence the need to move the
         // call to PrintStashedSupernovaDetails() in Star:EvolveOneTimestep())

src/BaseStar.h

@@ -191,7 +191,7 @@ class BaseStar {
 
 
     // setters
-            void                SetAngularMomentum(double p_AngularMomentum)                    { m_AngularMomentum = std::min(p_AngularMomentum, 0.0); }
+            void                SetAngularMomentum(double p_AngularMomentum)                    { m_AngularMomentum = std::max(p_AngularMomentum, 0.0); }
             void                SetInitialType(const STELLAR_TYPE p_InitialType)                { m_InitialStellarType = p_InitialType; }
             void                SetError(const ERROR p_Error)                                   { m_Error = p_Error; }
             void                SetObjectId(const OBJECT_ID p_ObjectId)                         { m_ObjectId = p_ObjectId; }

src/Star.cpp

@@ -176,6 +176,22 @@ bool Star::RevertState() {
     return result;
 }
 
+/*
+ *
+ * Resolve the loss of an envelope, including switching the stellar type.
+ *
+ * Keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency
+ * (envelope loss is rapid, so no time for angular momentum transport)
+ *
+ * void ResolveEnvelopeLossAndSwitch()
+ *
+ */
+void Star::ResolveEnvelopeLossAndSwitch() {
+    double omega = Omega();
+    (void)SwitchTo(m_Star->ResolveEnvelopeLoss(true));
+    SetOmega(omega);
+}
+
 
 /*
  * Apply mass changes if required, age the star one timestep, advance the simulation time, and update the

src/Star.h

@@ -242,7 +242,7 @@ class Star {
                                                    const double p_CompanionRadius,
                                                    const double p_CompanionEnvelope)                                { return m_Star->ResolveCommonEnvelopeAccretion(p_FinalMass, p_CompanionMass, p_CompanionRadius, p_CompanionEnvelope); } 
 
-    void            ResolveEnvelopeLossAndSwitch()                                                                  { (void)SwitchTo(m_Star->ResolveEnvelopeLoss(true));  SetOmega(Omega()); }                                        // keep core rotating at same angular frequency, no time for angular momentum transport on rapid envelope removal
+    void            ResolveEnvelopeLossAndSwitch();
 
     void            ResolveShellChange(const double p_AccretedMass)                                                 { m_Star->ResolveShellChange(p_AccretedMass); }
 

src/changelog.h

@@ -1385,7 +1385,9 @@
 //                                      - However, subsequently, only CHE stars are artificially kept in co-rotation with binary (ignoring angular momentum conservation) only if TIDES_PRESCRIPTION::NONE is used
 //                                      - Clean-up of BaseBinaryStar::ResolveMassChanges(): if m_Mass variable has already been updated (because ResolveEnvelopeLoss() has been called), no need to update attributes again
 //                                      - Associated code clean-up
+// 03.08.01   IM - Nov 19, 2024     - Defect repairs
+//                                      - Multiple rotation-related fixes to 03.08.01 (units, initialisation, min->max typo, retain omega on envelope loss)
 
-const std::string VERSION_STRING = "03.08.00";
+const std::string VERSION_STRING = "03.08.01";
 
 # endif // __changelog_h__