Implementing the Shikauchi+ (2024) core mass prescription for main sequence stars

[brcekadam]

This PR addresses the underestimation of convective core masses for mass-losing main sequence (MS) stars in COMPAS by implementing the Shikauchi et al. (2024) prescription (https://arxiv.org/abs/2409.00460). This is my first major PR and also my first C++ project, so I will greatly appreciate any feedback. I am happy to discuss and make changes as needed! This PR adds the following functionality:

1. Shikauchi+ (2024) core mass prescription:

• Implements equations (12) and (13) to model the evolution of convective core mass and central helium fraction for MS stars with MZAMS $\geq$ 10 Msol.
• Includes interpolation/extrapolation of coefficients across all metallicities.
• Introduces a new option --main-sequence-core-mass-prescription, which has the following modes: SHIKAUCHI (uses the new implementation), MANDEL (replaces the --retain-core-mass-during-caseA-mass-transfer option), NONE (no treatment for MS core mass).
• The core mass is stored internally in the m_MainSequenceCoreMass variable (replacing m_MinimumCoreMass), and it is currently not part of the output.
• Added a new global variable m_TotalMassLossRate, which is generally equal to -m_Mdot,
  but gets updated to the mass transfer rate (mass loss rate for donors and mass gain rate
  for accretors) during a mass transfer episode

2. Updated stellar tracks:

• Implements a new luminosity formula based on the Shikauchi+ (2024) prescription (equation A5), dependent on the central helium fraction and convective core mass. This formula is valid during core hydrogen burning (until HeliumAbundanceCore = 1 - Metallicity).
• Maintains Hurley's radius calculations, but updates the effective age at each timestep to ensure the main sequence duration aligns with the new prescription (age of 99% tMS is reached when central helium fraction of 1-Z is reached - this is when the MS hook begins).
• Smoothly connects the last point of core hydrogen burning (in radius and luminosity) with the first Hertzsprung gap point using linear interpolation to ensure the smoothness of stellar tracks.

3. Rejuvenation for MS accretors:

• Extends equations (12) and (13) to account for accretors by modeling rejuvenation (as new hydrogen is brought into the core). Here we assume a linear profile (in mass) between Yc and Y0, conservation of helium mass in the core before and after accretion, and that the accreted gas is pristine (i.e. has a helium fraction of InitialHeliumAbundance).

Things to keep in mind or improve:

• This prescription can be currently used only for stars with MZAMS >= 10 Msol due to limitations in the luminosity fits (this can be improved in the future)
• For MZAMS < 10 Msol, the MANDEL option is used as a fallback.
• I have also added a treatment for MS mergers that takes advantage of the new prescription (for now just using a simple assumption that the core mass of the merger product can be determined from equation (A3)).

Plots:

Using the new prescription leads to significantly higher core masses at TAMS for single stars, especially for the more massive stars that experience significant mass loss during their life on the MS.

HR diagrams showing the updated SSE tracks.

In this plot, I am varying the initial separation of the binary to test how different extent of case A mass transfer affects the core mass at TAMS. The retained core mass is higher when using the new prescription.

Here is a more detailed view of what happens to the stars during case A mass transfer when the new prescription is used. I am only showing what happens to the stars while the donor is on the main sequence.

[brcekadam]

Some things I wanted to ask:

1. I realized that I didn't use utils::Compare() anywhere, should all float comparisons be replaced with utils::Compare()?
2. I am initializing some relevant parameters in MS_gt_07.h. It probably makes more sense to do that in MainSequence.h, I just might need help with how to do that properly.

[jeffriley]

@brcekadam I'll review in detail a bit later, but for now:

I realized that I didn't use utils::Compare() anywhere, should all float comparisons be replaced with utils::Compare()?

Mostly. This is why we use utils::Compare():

Compiling and running this code:

#include <iostream>
#include <iomanip>

int main(int argc, char* argv[]) {

    double a = 0.1;
    double b = 2.1;
    double c = 31.7;
    double d = 18.53;
    double e = 10.9;

    double v = a;
    
    std::cout << "@1: " << (v == a ? "v == a" : "v != a") << "\n";
    
    v *= b;
    v += c;
    v /= d;
    v -= e;

    v += e;
    v *= d;
    v -= c;
    v /= b;
    
    std::cout << "@2: " << (v == a ? "v == a" : "v != a") << "\n";

    std::cout << "a = " << a << ", v = " << v << "\n"; 
    std::cout << std::fixed << std::setprecision(15) << "a = " << a << ", v = " << v << "\n"; 

    return 0;
}

gives:

@1: v == a 
@2: v != a
a = 0.1, v = 0.1
a = 0.100000000000000, v = 0.099999999999995

With every manipulation we lose precision (or maybe accuracy...) - digits roll off the end of the number. And, as it happens, because of the way C/C++ stores floating-point numbers (IEEE standard), not all real values can be represented exactly.

utils::Compare(), when it is implemented properly - we haven't settled on tolerances yet, will compare with tolerance in an attempt to avoid problems like that shown above.

However, there are some exceptions throughout the code (should all be commented). The main reason for not using utils::Compare() is when you just want to know the sign of a value (to prevent (say) sqrt of a negative number) - then you would use e..g if (x >= 0.0) y = sqrt(x) rather than utils::Compare().

I am initializing some relevant parameters in MS_gt_07.h. It probably makes more sense to do that in MainSequence.h,
I just might need help with how to do that properly.

I'll take a look at that when I review later

[ilyamandel]

Thank you for the very nice work, @brcekadam !

Some of the requests below may overlap with @jeffriley 's.

For now, I am just reviewing the class with most of the changes, MainSequence.cpp.

m_HeliumAbundanceCore -- does this hold a meaningful value only if the Shikauchi prescription is used? If so, perhaps reflect that in the name so another developer doesn't misuse it?

I recommend defining the 10.0 Msun threshold in constants.h -- this way, if someone changes it, they'll automatically consistently change it everywhere it's used.

Around line 350:
if ((p_HeliumAbundanceCore == 1.0 - m_Metallicity) && (m_TotalMassLossRate == -m_Mdot))
You probably want utils.Compare() >= 0.0 on the first test; perhaps add a comment that, once you get to the hook, the He abundance fraction will be fixed to 1-metallicity and the core is not evolving-- otherwise, not clear why the luminosity being computed from Eq. (A5) when this function is called, say, half-way through the hook can be used as the luminosity at the start of the hook.
For the second one, utils.Compare == 0.0 ?
(see also line 538)

--
CalculateCoreMassMandel()
There's a function called TAMSCoreMass(); use that instead of cloning, so you just need a one-liner:

return std::max(m_MainSequenceCoreMass, CalculateTauOnPhase() * TAMSCoreMass() );

double totalMass             = m_Mass;
double centralHeliumFraction = m_HeliumAbundanceCore;
double mixingCoreMass        = m_MainSequenceCoreMass;
double lnMixingCoreMass      = std::log(mixingCoreMass);
double heliumFractionOut     = m_HeliumAbundanceCoreOut;

Other than the line which computes a log (which you obviously don't want to call more than once), I really don't see the point of renaming variables -- just makes the code longer and harder to read... (and even the ln you only use once, so no point in defining a new variable)

Line 815:
// Eq (A4)
double beta = 1.0 - FMIX_COEFFICIENTS[1] * totalMass / (FMIX_COEFFICIENTS[2] * fmix) * std::exp(-totalMass / FMIX_COEFFICIENTS[2]);
Eq. (A4) refers to M_ZAMS; why does this use m_Mass?

Line 831:
I don't understand what you are doing with the ODE solver. You aren't actually changing the value of \dot{Y} -- your right-hand-side of dxdt[0] is constant -- so why evolve it if it's a constant? If you do want to evolve it, luminosity and mass are changing as functions of time, too -- e.g., the mass is not mixingCoreMass, its exp(x[1]) -- but do you need to take small internal step sizes here since we are just doing difference equations elsewhere with the step size p_Dt? Except, perhaps, fast case A MT -- but then, the right-hand-side shouldn't be constant, either...

You don't need to muck with time units -- Q_CNO cannot be LSOL/MSOL (it should be energy per mass, contrary to what constants.h claims), so just define it in units of LSOL * Myr / MSOL, and then you can keep timesteps in their code units.

Oh, and don't take log10(luminosity) just to take PPOW(10.0,logL) a moment later -- it's unnecessary computational cost. :)

Line 855:
I don't understand what you are integrating. As far as I can tell, the right hand side is a constant:
dxdt[0] = (heliumFractionOut - m_InitialHeliumAbundance) / (mixingCoreMass - m_InitialMainSequenceCoreMass) * (deltaCoreMass / currentTimestepInYrs);
If dx/dt = k, where k is a constant, there is no need to integrate numerically :-) -- the solution is x(t) = x(0) + k t.
It's also possible I missed something -- adding comments would definitely help to clarify what you are doing in places like this, especially if there are no published equations to refer to.

[I won't comment further on rejuvenation because I couldn't quite follow what you were doing here; can I suggest that you write up some notes on the equations you are trying to solve in latex and attach them to this PR as a PDF? You'll reuse them later for your paper...]

Line 929 & 932: another example of unnecessary variable declaration
mixingCoreMass = std::get<0>(ShikauchiSolution);
m_MainSequenceCoreMass = mixingCoreMass; // Update the core mass

Around lines 927--935: the if and else if clearly don't cover all of the cases, and I am a bit worried that we are doing nothing when neither set of conditions are satisfied...

[ilyamandel]

@reinhold-willcox 's comment today about implementing wind mass accretion is one of the reasons why
m_TotalMassLossRate == -m_Mdot
in MainSequence::CalculateLuminosityShikauchi()
is a dangerous way to check whether mass transfer is ongoing: the mass may be changing because the star is accreting winds, even in the absence of mass transfer...

[ilyamandel]

MainSequence::CalculateLuminosityShikauchi()
if ((p_HeliumAbundanceCore == 1.0 - m_Metallicity)
I think it would be safer to do if(utils::Compare(p_Age, 0.99 * tMS) >= 0)

I don't quite understand why MainSequence::CalculateMainSequenceCoreMassShikauchi() returns the new central He fraction and core mass, but updates m_HeliumAbundanceCoreOut directly. That seems inconsistent. Why not just update all of the variables directly (and make the function return void)?

In the same function, you don't really need an ODE solver for Eq. (7) in your notes, especially since you are treating Eq. (5) as exact.

Same function: if a star accretes a significant amount of mass and go from M_c<M_c,0 to M_c>M_c0 in one time step, the derivative term should be used for all the mass accreted up to M_c,0 but not beyond. (See Slack discussion and @ryosuke-hirai 's comment on this.)

void MainSequence::UpdateMainSequenceCoreMass()
ShikauchiSolution = CalculateMainSequenceCoreMassShikauchi(p_Dt);
mixingCoreMass = std::get<0>(ShikauchiSolution);
centralHeliumFraction = std::get<1>(ShikauchiSolution);
m_HeliumAbundanceCore = std::min(centralHeliumFraction, 1.0 - m_Metallicity); // Update the core helium abundance
m_MainSequenceCoreMass = mixingCoreMass; // Update the core mass
As mentioned previously, just use the stored values instead of passing them around. Also, CalculateMainSequenceCoreMassShikauchi() should internally cap the central helium fraction to 1.0-m_Metallicity.

            p_TotalMassLossRate != m_TotalMassLossRate

That's a very non-intuitive comparison. :)

[brcekadam]

Thank you so much for your comments, @jeffriley and @ilyamandel! I have now fixed most things, but here are some comment and answers to your questions:

--main-sequence-core-mass-prescription MANDEL has the exact same functionality as --retain-core-mass-during-caseA-mass-transfer, and the detailed output should be the same. I put the --retain-core-mass-during-caseA-mass-transfer option back as a deprecated option.

I have restructured the conditions in MainSequence::UpdateMainSequenceCoreMass(), but not doing anything for some cases is intentional. I only want to do the core mass calculations if time advances, star is not in the MS hook (the prescription does not apply for those), and I want to prevent the function from being called as part of SSE if it was already called as part of BSE for that time step (more on that below).

m_HeliumAbundanceCore does have a meaningful value without Shikauchi prescription being used. A linear model is implemented as a default, in which core helium abundance linearly increases with m_Tau.

I have defined the threshold for Shikauchi prescription in constants.h and changed it to 15 Msun for now.

Regarding the m_Mass used in eq. (A4) instead of m_MZAMS, if you look at equation (12) in the paper, current mass is used as the argument in function Beta. I checked with @ryosuke-hirai before, and it should be the current mass that’s used.

I've changed what is being integrated in the ODE solver. I was indeed integrating a constant term, which works okay for small mass loss (e.g. through winds). But for thermal timescale mass loss, we definitely need to take time steps smaller than the COMPAS time step to fully resolve the core mass evolution.

This comparison p_TotalMassLossRate != m_TotalMassLossRate is something that I don’t really know how to avoid. I changed the name of the function argument to p_MassLossRate to make it a little less confusing at least. My issue is that when UpdateMainSequenceCoreMass() gets called during binary evaluation (when there is ongoing mass transfer), I want to prevent it from being executed again in SSE when the individual stars are aged, since this would cause the calculation to be executed twice for a single time step. This is the only reason why I even introduced the m_TotalMassLossRate variable. It is generally equal to -m_Mdot (it gets updated whenever BaseStar::CalculateMassLossRate() is called), but gets updated to the mass transfer rate during a mass transfer episode, and it stays as that value when UpdateMainSequenceCoreMass() gets called as part of SSE, so p_MassLossRate != m_TotalMassLossRate. This approach works, since we are not considering wind mass loss during mass transfer, but I can also see how this could cause problems in the future. I couldn’t really figure out a better way to do it, but I am open to suggestions if you have any ideas :-)

I am also attaching a PDF document describing the rejuvenation procedure.
Rejuvenation_in_main_sequence_accretors (3).pdf

[jeffriley]

@brcekadam I'm working on this - I will do my best to have it done before the end of this week. There are some changes I'd like to see, and I'll probably push those myself (it'd take me longer to describe them than actually make the changes). The initialisation of the constants bothers me a bit - I'm trying to decide the best place for that. I might suggest moving them to BaseStar and only do the initialisation if the option is enabled - but I'm not convinced I'm comfortable with that either. The problem is that because they are in MS_gt_07 they rely on the SHIKAUCHI threshold always being above 0.07 MSun - which I'm sure it always will be, but me (and you) being sure is not the same as it being that way in the code... Make sense? So, let me think about that and I'll get this done in the next day or two...

[ilyamandel]

@jeffriley -- I'd be OK with a big fat comment warning or an assert statement that SHIKAUCHI_LOWER_MASS_LIMIT must be above 0.7 Msun. ;-)

[ilyamandel]

Very nice, @brcekadam ! One very minor request, one query, and one idea to consider. I am listing this as an approval, so you don't need to wait for another review from me. But, of course, I'd be interested in hearing your thoughts -- and please do wait for Jeff's comments.

• the comment on line 761 of MainSequence.cpp is missing one of the function arguments:

• void CalculateMainSequenceCoreMassShikauchi(const double p_Dt)

• I didn't know you could declare variables to have type auto in C++ (e.g., around lines 776-780 of MainSequence.cpp). Cool, thanks for teaching me! For my education, in this context, what's the benefit of doing that rather than declaring them to be double? [This might be a question for both you and @jeffriley :) ]

On the p_TotalMassLossRate issue:
Actually, I think the clean/safe thing to do is to call the function twice. :-)
Once for RLOF, but only taking into account the RLOF mass loss/gain; then again for winds, only taking into account the winds mass loss (or gain, if we ever implement it). Of course, only the second one would ever get called in SSE...

[jeffriley]

@ilyamandel

I'd be OK with a big fat comment warning or an assert statement that SHIKAUCHI_LOWER_MASS_LIMIT must be
above 0.7 Msun. ;-)

That might end up being the solution, but I'm not yet convinced it is the best one. I'm struggling with functions that rely on the star being a MS_gt_07 star not being in that class (they're in the MS class). I get why @brcekadam has done it the way he has, I'm just not sure it's the best way. But fear not, I will be pragmatic...

I didn't know you could declare variables to have type auto in C++

We use auto fairly extensively throughout the code (grep for 'auto'...). I use it mostly for iterators when I don't want to make a mistake with the type of iterator - so I let the compiler deduce it :-)

@brcekadam uses it here to type the return value of nested functions. There is no real reason to use auto here - all those functions will return doubles, so the return values should really be declared as doubles (if it should be a double, type it as a double - if it could be something else, maybe the code needs to change...). That's probably true for the ode too, but I haven't actually looked at that yet.

EDIT: actually, I might take that back :-) I think auto might be correct there (because it is a nested function) - I'll take a closer look (I think it depends where it is called...)

EDIT: Ah, yes - the nested functions are lambda functions. They could probably be rewritten as properly typed and named nested functions, but lambda functions are ok here. So, I take it back :-)

[brcekadam]

Thanks @ilyamandel!

They are indeed lambda functions. Alpha, Beta, and Delta are functions of the core mass, total mass, and central helium fraction, so I thought this would be a clean way of defining them as nested functions. They are not needed anywhere else outside CalculateMainSequenceCoreMassShikauchi().

Regarding the p_TotalMassLossRate issue... I am not sure if we could do that. The problem is that equation (12) has two terms: one for the natural decay of the core and second term reflecting the effect of mass loss. If I do call the function twice in one time step (once in BSE and once in SSE), I would need to set the natural decay term to zero either in BSE or SSE to avoid unphysical ageing of the core. The natural decay term is important in BSE for the rejuvenation calculations to work properly and I also can't set it to zero in SSE, so this is a bit tricky. Ideally, we would want to call the function only once per time step with all the mass loss/gain accounted for.

[jeffriley]

@brcekadam

I have some changes that I can't push at the moment due to a github outage - I'll push them as soon as I can. Here are my comemnst (associated with the changes I will push):

What does CORE_MASS_PRESCRIPTION::NONE really mean? It looks like it really means CORE_MASS_PRESCRIPTION::ZERO. Should we rename it to be consistent with other prescriptions that have "ZERO" options (rather than "NONE")?
(I think "NONE" for a prescription means "we don't do anything for this case - it's a no-op", whereas "ZERO" means "we set the value to 0.0"). Hmmm - here that might be ambiguous - no core mass treatment probably means core mass = 0.0... Your call :-)

Wrt MainSequence::CalculateMainSequenceCoreMassShikauchi(), most "Calculate" functions just calculate and return (there are a couple of exceptions that we should address one day...). I changed this one to do the same.

I removed the class variable m_HeliumAbundanceCoreOut, and replaced it with a static local variable in MainSequence::CalculateMainSequenceCoreMassShikauchi() (heliumAbundanceCoreOut ).

Can we do the same with m_InitialMainSequenceCoreMass? As far as I can tell we only use it in MainSequence::CalculateMainSequenceCoreMassShikauchi(), except that we reset it in MainSequence::UpdateAfterMerger().

I feel like we could do it if we move this clause in MainSequence::UpdateAfterMerger():

if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) {
    m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass);           // update initial mixing core mass
    m_MainSequenceCoreMass        = m_InitialMainSequenceCoreMass;                          // update core mass
}

into MainSequence::CalculateMainSequenceCoreMassShikauchi(), and call MainSequence::CalculateMainSequenceCoreMassShikauchi() from MainSequence::UpdateAfterMerger().

In MainSequence::CalculateMainSequenceCoreMassShikauchi() we would declare a new variable (to replace m_InitialMainSequenceCoreMass), initialMainSequenceCoreMass as static, and key off the fact that m_MainSequenceCoreMass was just set = 0.0 in MainSequence::UpdateAfterMerger(), then something like:

static double heliumAbundanceCoreOut{ m_InitialHeliumAbundance };                                                                                           // Helium abundance just outside the core - initially m_InitialHeliumAbundance
static double initialMainSequenceCoreMass{ CalculateInitialMainSequenceCoreMass(m_MZAMS) };
if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0) && utils::Compare(m_MainSequenceCoreMass, 0.0) == 0) {
    initialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(m_Mass);           // update initial mixing core mass (m_Mass just set in MainSequence::UpdateAfterMerger())
return std::tuple<double, double> (initialMainSequenceCoreMass, heliumAbundanceCoreOut);
}

Would that work?

If we can remove the class variable m_InitialMainSequenceCoreMass (and make it local to MainSequence::CalculateMainSequenceCoreMassShikauchi()), then MainSequence::CalculateMainSequenceCoreMassShikauchi() could be const.

MainSequence::CalculateInitialMainSequenceCoreMass() description, and description of parameter and return value, is not strictly correct - the parameter is not always the ZAMS value (sometimes it is post-MS merger)

MainSequence::CalculateConvectiveEnvelopeMass() return type and return value in description not same as code (returns a tuple)

I am not completely happy with the initialisations done in MS_gt_07.h, but for now I don't have a better solution (maybe we should just put it in BaseStar...). I'm uneasy because the code as it is relies on MS:_gt_07::Initialise() being executed - but when we implement starting evolution at stellar types after ZAMS, MS:_gt_07::Initialise() won't be executed. Will that matter? It doesn't right now, but it just makes me uneasy that the dependency is there and may be forgotten as time passes...

[jeffriley]

@brcekadam @ilyamandel I just pushed my changes

[ilyamandel]

@brcekadam -- yup, I see your point. And you can't expect to be able to access variables like m_MassLossRateInRLOF or functions like IsRLOF() because you don't necessarily have a binary star. OK, maybe this is the best we can do, though it feels like it may prove difficult to maintain in the future...

[brcekadam]

Thank you, @jeffriley!

You're right, I'll change the name to CORE_MASS_PRESCRIPTION::ZERO. It will be more consistent with other prescriptions, and we are setting the value to zero when this prescription is used, so it does make more sense :-)

I actually wasn't familiar with static variables! I was doing a bit of testing and if I understand it correctly, heliumAbundanceCoreOut gets updated to a new value every time MainSequence::CalculateMainSequenceCoreMassShikauchi() is called, right? I noticed that what ends up happening is that if we have a binary and Star1 is aged, MainSequence::CalculateMainSequenceCoreMassShikauchi() is called, heliumAbundanceCoreOut is updated to a new value, and it is this value that gets passed to Star2 as the initial value when it's being aged. After being updated, it gets passed to Star1 as the initial value, and so on. So this one heliumAbundanceCoreOut is being passed around between the two stars, which could cause some strange behaviour. We would need each star to keep track of its own heliumAbundanceCoreOut and initial core mass, is there a way to do that with static variables?

Regarding the initialisation... I'll keep this in mind when starting evolution at different stellar types is implemented. I'll have to change a few things anyway, since we don't necessarily know what the core mass at TAMS is without fully evolving the star through the MS. Maybe I could do a fit of the data that's shown in the first plot (core mass at TAMS vs. initial mass), so that we know what the core mass at TAMS is for a given initial mass, and then that core mass value will be initialised post-MS. I'll be happy to work on this in the future. One thing I was also thinking about it is how the initialisation in MS_gt_07.h works in the context of CH stars. When CH stars switch to MS after spinning down, is MS_gt_07::Initialise() called? We probably don't want to reset the core mass to the initial core mass in that case. Would initialising in MainSequence.h instead of MS_gt_07.h help with this?

[jeffriley]

I noticed that what ends up happening is that if we have a binary and Star1 is aged ...

Ah, yes - my mistake. In c++, a static local variable belongs to the class, not an instance of the class - so the behaviour you are seeing is expected (all instances of a class share the same static local variable) (I had forgotten that - I have always thought the implementation is not right - I think it makes more sense for the static local variable to belong to the instance, not the class, but apparently the c++ standard developers don't agree with me...).

We can get around it though :-) In the code below I declare a static local variable that is actually a map, where the key is the instance id and the value is the value we want to store (in the code below, just an integer). Accessing the map variable using the instance id (this) accesses the correct element of the map. It's kind-of a kludge, but not really - we're just using a valid feature of the language to our advantage :-) Since binaries only have two stars the overhead is minimal, and it lets us use static local variables (in the way I think they should be :-))

#include <iostream>
#include <map>

class thisClass {
    public:
    void incrementAndPrint(const std::string str) {
        static std::map<thisClass*, int> v;
        ++v[this];
        std::cout << str << "::v = " << v[this]<< '\n';
    }
};

int main(int argc, char* argv[]) {

    thisClass a1, a2;

    a1.incrementAndPrint("a1");
    a1.incrementAndPrint("a1");
    a2.incrementAndPrint("a2");
    a1.incrementAndPrint("a1");
    a2.incrementAndPrint("a2");
    a1.incrementAndPrint("a1");
    
    return 0;
}

EDIT: the code above produces the following output

a1::v = 1
a1::v = 2
a2::v = 1
a1::v = 3
a2::v = 2
a1::v = 4

end edit

(Aside: the way static local variables are implemented in c++ does open up the possibility that constituents of a binary star could actually share data - that could be interesting (even just to share signals/messages)).

When CH stars switch to MS after spinning down ...

Unless @ilyamandel or @veome22 changed this (and I confess I don't remember), CH stars don't spin down - or at least they don't change stellar type if they do.

EDIT:

Since binaries only have two stars the overhead is minimal

Clones could affect that, but still shouldn't impose too much overhead - we can keep an eye on that.

[jeffriley]

Hmmm. Actually... I suspect the map in the code above might be used for every star/binary in the run... That could impose a bit more overhead than I thought...

Maybe going back to using a class variable is better. We could test performance I suppose, but I think a class variable is safer. Bummer...

[brcekadam]

@jeffriley @ilyamandel

I have just pushed the final changes! I renamed the NONE prescription to ZERO, and brought back the class variable m_HeliumAbundanceCoreOut. The issue that came up at the end of our code call yesterday was not an issue at all actually. I wasn't running the latest version of COMPAS, it was saying that I was on version 3.11. because I set it to be that in changelog.h, but I never pulled the latest changes. So all is good, now I am getting what you were getting :-)

[jeffriley]

All good - thanks @brcekadam !

[ilyamandel]

Unless @ilyamandel or @veome22 changed this (and I confess I don't remember), CH stars don't spin down - or at least they don't change stellar type if they do.

I can't find anywhere in the code where we distinguish between CHE_MODE::PESSIMISTIC and CHE_MODE::OPTIMISTIC. We clearly should, but you are right, we seem to ignore this setting! That's a bug, right, @jeffriley?

[jeffriley]

Yes, it would seem it is a bug... I'm sure that functionality was once in the code...

…

On Fri, 17 Jan 2025, 5:24 pm Ilya Mandel, ***@***.***> wrote: Unless @ilyamandel <https://github.com/ilyamandel> or @veome22 <https://github.com/veome22> changed this (and I confess I don't remember), CH stars don't spin down - or at least they don't change stellar type if they do. I can't find anywhere in the code where we distinguish between CHE_MODE::PESSIMISTIC and CHE_MODE::OPTIMISTIC. We clearly should, but you are write, we seem to ignore this setting! That's a bug, right, @jeffriley <https://github.com/jeffriley>? — Reply to this email directly, view it on GitHub <#1292 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGM5XNSGHZCXDA5PLZ3PMU32LCOY3AVCNFSM6AAAAABST6UFESVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKOJXGUZDGNJQGQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[ilyamandel]

@jeffriley : I created a new issue, #1323 .