Discontinuous radius evolution

ilyamandel commented 1 year ago

We still have an issue with discontinuous radial evolution in COMPAS.

Here's a trivial example, from

./COMPAS -n 1 --mode SSE --initial-mass 50 --detailed-output

Note that this issue does not come up in Hurley's SSE/BSE.

I believe this may be due to differences between the SSE paper and its implementation.

For example, around line 190–200 of hrdiag.f:

rmin = rminf(mass)

This uses mass (i.e., M_0 in the language of the paper) to calculate rmin; meanwhile, the paper states that M_0 is used for timescales and luminosities, but M (mt in the code) for radii.

ry = ragbf(mt,lums(4),zpars(2))
rx = MIN(rmin,ry)
…
tau2 = tblf(mass,zpars(2),zpars(3))
if(tau2.lt.tiny) rx = ry

This means that, for massive stars, rmin isn’t used anyway — instead, it’s the AGB luminosity that’s used directly for the luminosity at the end of the Hertzsprung gap. And that fundamentally explains how luminosities are consistent between HG and CHeB stars even though, according to the formulae in the paper (as I read them), one uses a function of mass and the other of mass and luminosity.

In the end, it makes sense: the HG is just a smooth connection from the MS that precede it to the giant phases that follows. We need to make sure that our code does the same.

reinhold-willcox commented 1 year ago

Good catch! Have you talked to Jarrod about why there is a discrepency?

ilyamandel commented 1 year ago

@reinhold-willcox -- yes, the discrepancy is there because the SSE paper doesn't quite cover match the SSE code in some cases; the trouble for me at the moment is parsing the SSE code (done for HG stars, but still struggling for CHeB stars, see Slack).

jeffriley commented 6 months ago

@ilyamandel

I have made a couple of changes that I think help fix the problem - I can get COMPAS to look like Hurley sse for a single star (haven't tried others) with mass loss disabled. Happy to try with mass loss enabled if you tell me how to configure sse for mass loss and tell me which option(s) in COMPAS match.

Branch on my fork - hopefully you have access (let me know if not). Look for "~~ILYA~~" in the code to see the changes made. Explanatory comments there.

https://github.com/jeffriley/COMPAS-1/tree/hurley-test

I modified Hurley sse to print a "detailed output file":

I put this at the top of evolve1.f (just after the COOMON blocks):

OPEN(73,file='detailed.dat',status='unknown')

and this just after the call to hrdiag at line 169:

          write(73,7373) aj, dtm, tphys, kw, mt, r, lum
7373  format(f25.15, ',', f25.15, ',', f25.15, ',', i2, ',',
        &       f25.15, ',', f25.15, ',', f25.15)

(For the files below, you'll need to remove the .txt suffix - that's there to convince github to allow the uploads - except for the timesteps file - that one can stay)

I used this evolve.in file:

evolve.in.txt

and the run resulted in this "detailed output file":

detailed.dat.txt

I took the timesteps from the sse detailed output and produced this timesteps file for COMPAS:

sse_timesteps_1.txt

and ran COMPAS with this command using the sse timesteps:

./compas --mode sse -n 1 --random-seed 0 --initial-mass 50 --metallicity 0.02 --use-mass-loss false --timesteps-filename sse_timesteps_1.txt --detailed-output --logfile-type csv

and produced this CSV detailed output file (don't forget to filter by record type, and only look at record type 4):

SSE_Detailed_Output_0_sse_timesteps.csv

and also ran COMPAS with this command and let COMPAS determine the timesteps:

./compas --mode sse -n 1 --random-seed 0 --initial-mass 50 --metallicity 0.02 --use-mass-loss false --detailed-output --logfile-typ

e csv

and produced this CSV detailed output file:

SSE_Detailed_Output_0_compas_timesteps.csv

The COMPAS detailed output file where COMPAS used the sse timesteps matches the sse output exactly (except that COMPAS doesn't evolve the black hole). The detaield output file where COMPAS determined the timesteps ends up at the same place, but we have a lot more timesteps on the CHeB phase.

Let me know what you think. More test required, but I think this is a step forward.

ilyamandel commented 6 months ago

Thank you, @jeffriley , that looks very promising!

I think it may be tricky to get SSE and COMPAS mass loss prescriptions to match. However, we can test for whether single-star radial evolution in COMPAS, after allowing for mass loss, would be smooth. The original post in this thread had a simple test case and showed a plot with non-smooth radial evolution. How does your fixed code perform in this case?

./COMPAS -n 1 --mode SSE --initial-mass 50 --detailed-output

jeffriley commented 6 months ago

I have done some more work on this to see if we can get COMPAS to match Hurley's sse with mass loss enabled.

I have identified a number of places in the Hurley sse code where the code differs from the description of the algorithm in the Hurley et al. 2000 paper. These are:

The helium rate constant defined in Hurley et al. 2000, eq 68 as:

AHe = (EXHe)−1 = 7.66 × 10−5M⊙ L⊙−1Myr−1

The value used in the code (star.f, lines 76 & 318) as 8.0E-05. COMPAS uses 7.66E-05.
The value of the solar temperature in kelvin used in the Hurley sse code (evolv1.f line 266) is calculated as:

1000.0 * (1130.0 ^ 0.25) = 5797.885. COMPAS uses the value 5778.0.
The radius calculation on the HG in the Hurley sse code (hrdiag.f lines 92 & 188-203) differs significantly from Hurley et al. 2000, eq 27, and the COMPAS code (which matches Hurley et al. 2000, eq 27, see HG::CalculateRadiusOnPhase()) Also, hrdiag.f lines 188-203 indicate that we should use M0 to calculate the radius on the HG, not M (COMPAS uses M, see HG::CalculateRadiusOnPhase() declaration in HG.h)
The luminosity calculation on the CHeB phase in the Hurley sse code includes a check for high-mass stars (hrdiag.f lines 287-302) that is not present in the COMPAS code (see CHeB::CalculateLuminosityOnPhase()). Luminosity on the CHeB phase is calculated using Hurley et al. 2000, eqs 61, 62 & 63, but I don't see a check for high-mass stars anywhere in those equations or surrounding discussion. There does seem to be a typo in that discussion:

"Hence, τy = 1 for LM and IMstars and τy = τbl for IMstars".
Hurley et al. 2000, just after eq 105 has the following text:

"If the star is on the HG or GB then we have, for M < MHeF ..."

This is a typo - it should read "M > MHeF". Additionally, the Hurley ss code uses mass0 for the check, not mass. This can be verified by inspecting the Hurley sse code, hrdiag.f, line 702.

From hrdiag.f:

*
- Calculate the core radius and the luminosity and radius of the
- remnant that the star will become.
- ```
 tau = 0.d0
 if(kw.le.1.or.kw.eq.7)then
    rc = 0.d0
 elseif(kw.le.3)then
    if(mass.gt.zpars(2))then            <----- line 702
       lx = lzhef(mc)
       rx = rzhef(mc)
       rc = rx
    else
       if(wdflag.eq.0)then
          lx = 635.d0*mc*zpars(14)/((ahe*0.1d0)**1.4d0)
       elseif(wdflag.ge.1)then
          lx = 300.d0*mc*zpars(14)/((ahe*0.1d0)**1.18d0)
       endif
       rx = 0.0115d0*SQRT(MAX(1.48204d-06,
&           (mch/mc)**(2.d0/3.d0)-(mc/mch)**(2.d0/3.d0)))
       rc = 5.d0*rx
    endif
 elseif(kw.eq.4)then
```
The affects GiantBranch::CalculateRemnantRadius() and GiantBranch::CalculateRemnantLuminosity() in COMPAS, which both use "<" instead of ">"

I have made the following changes to the COMPAS code to address the issues noted above:

No changes to COMPAS (see note below re tests).
No changes to COMPAS (see note below re tests).
HG::CalculateRadiusOnPhase() rewritten to match Hurley sse code - see hrdiag.f lines 92, 188-203. HG::CalculateRadiusOnPhase() in HG.h changed to use m_Mass0 instead of m_Mass (per hrdiag.f lines 188-203).
Added check for high-mass star to CHeB::CalculateLuminosityOnPhase() - see Hurley sse hrdiag.f line 287.
GiantBranch::CalculateRemnantRadius() and GiantBranch::CalculateRemnantLuminosity() both changed to use ">" instead of "<".

I have made the following changes, either to support the changes noted above (e.g. new comstants in constants.h), or to fix issues noticed whilst making those changes:

I changed the value of the constant LOG10_ZSOL in constants.h from -1.69897 to -1.698970004336019. I made this change because the constant LOG10_ZSOL is used heavily in calculating the a and b coefficients (see Hurley et al. 2000, Appendix, page 24). Increasing the precision of LOG10_ZSOL results in subtle changes to the a and b coefficients, and in-turn subtle changes to the age of stars and timescales calculated using the coefficients.
Added EDDINGTON_PARAMETER_FACTOR constant in constants.h. The value is HE_RATE_CONSTANT * 0.325. As noted above, HE_RATE_CONSTANT is defined as 7.66E-05 in constants.h. I added this constant because the value 7.66E-05 * 0.325 is used in several places in the COMPAS winds code in BaseStar.cpp. If the use of 7.66E-05 was just a coincidence there and doesn't actually relate to the HE_RATE_CONSTANT constant, then we should change the definition of the EDDINGTON_PARAMETER_FACTOR constant in constants.h so we know what the value refers to.
Added HIGH_MASS_THRESHOLD constant in constants.h. The value is 12.0. This supports checks in the code for high-mass stars (some existing, some new as noted above).
A added the "virtual" keyword to the declarations of CalculateMassLossRate() and CalculateMassLossRateWolfRayetShenar2019() in BaseStar.cpp Both these fuctions have counterparts declared in derived classes. I haven't checked what the impact of this might have been - probably minimal (though @jmerritt1 and @SimonStevenson may want to check for CalculateMassLossRateWolfRayetShenar2019() since it was a recent addition for the newly updated/added winds prescriptions. CalculateMassLossRate() has counterparts in the NS and BH classes - I'm not sure when they were added (might have been me a while ago...), but it looks like they were added as insurance, and the impact there is probably minimal. The consequence of not including the keyword virtual on the base class declaration is that any counterpart functions in derived classes will be statically bound at compile time, thus possibly reducing runtime polymorphism - and perhaps not calling the functions as intended.

The changes noted above will be in COMPAS v02.44.01 if the open PR is merged before any other PRs are merged.

I ran some tests to verify the changes. The goal was to have COMPAS evolve single stars as close as possible to the Hurley sse evolution of single stars. I used COMPAS to evolve a number of single stars, each using the following COMPAS options:

`--mode sse`
`-n 1`
`--random-seed 0`
`--metallicity 0.02` 
`--mass-loss-prescription HURLEY`

The initial mass of each star was varied via the --initial-mass option. Stars were evolved with the following initial mass (Msol): 1, 3, 8, 10, 15, 25, 50, 80, 100, 120, and 150.

To compare agains Hurley sse evolution, I evolved a matching set of stars (i.e. the same initial masses) using the Hurley sse code, with the following changes:

I changed the helium rate constant used in the Hurley sse code to 7.66E-05 to match the value used in COMPAS.
I changed the solar temperature used in the Hurley sse code to 5778K to match the value used in COMPAS.
I changed the hewind parameter in the evole.in file from 0.5 to 1.0. hewind is used to calculate the small envelope parameter mu described in Hurley et al. 2000, eqs 97 & 98. COMPAS uses mu = 0.0 (see HeMS::CalculateMassLossRateHurley()) as indicated in Hurley et al. 2000 in the discussion after eq 106 (towards the bottom of the left column). As distributed, the evolve.in file used by the Hurley sse code uses a value of 0.5 for the hewind parameter, which would mean the value of mu would be 0.5. Comments in the Hurley sse code (see sse.f) indicate that the default hewind value should be 1.0, and since that agrees with comments in Hurley et al. 2000, and with the value used in COMPAS, I set hewind to 1.0 in the evolve.in file.

The results of the tests are very promising. Apart from very minor differences, probably due to differences in timestep calculations, and perhaps one issue that we may have to follow-up (see the discussion of the 8 Msol star test below), with the changes to COMPAS described above, evolution of single stars in COMPAS, with or without mass loss, matches sse evolution very closely.

The test run for a star with an initial mass of 8 Msol showed that while the sse code evolved the star through the MS, HG, FGB, CHeB, EAGB, and TPAGB before terminating evolution with a ONeWD, COMPAS missed the TPAGB phase entirely and evolved directly from an EAGB star to a ONeWD. This may have been because of small differences in the star's attributes that determine what phase it evolves to after the EAGB phase, but will need more investigation. It should be noted that this behaviour is evident in the current version of COMPAS (v02.44.00), so is not a regression caused by the changes described above (so maybe we just haven't noticed it before...).

I ran all COMPAS runs with --timestep-multiplier = 0.1. I did this because an initial run for a 50 Msol star indicated the COMPAS timestep calculations result in poor resolution on some phases (and too much resolution on others) - it seems we take too few, and large, timesteps on some phases, and probably too many, and too small (i.e. we waste compute cycles) on others.

This is the HR diagram for a 50 Msol star for sse and COMPAS at the default timestep multiplier (i.e. 1.0)

HRD

and with a timestep multiplier of 0.1:

HRD

I have produced some comparison plots for each of the tests run. The plots show the HR diagram, the evolution of mass, mass0, and core mass over time, and the evolution of radius over time. I have shown individual plots for both sse and COMPAS, as well as overlaid plots (the individual plots allow closer inspection of areas of small differences that may be obfuscated in the overlay plots). Each pdf file below contains the plots for a single test (different initial masses). Full resolution plots, and the datafiles used (detailed output files for both COMPAS and sse - filenames describe the test: the sse filenames include the initial mass value, and the COMPAS filenames include the initial mass value and the timestep-multiplier value), can be found at

https://drive.google.com/drive/folders/1GlBi-cbVUY1Rnlcpe3b_WE33NITg2epG?usp=drive_link

1 Msol: plots.pdf 3 Msol: plots.pdf 8 Msol: plots.pdf 10 Msol: plots.pdf 15: Msol: plots.pdf 25 Msol: plots.pdf 50 Msol: plots.pdf 80 Msol: plots.pdf 100 Msol: plots.pdf

The maximum initial mass allowed by sse is 100 Msol, so the following two sets of plots show results for COMPAS only

120: Msol plots.pdf 150 Msol: plots.pdf

jeffriley commented 6 months ago

See PR #1102

ilyamandel commented 6 months ago

Fixed in #1102 , closing.

TeamCOMPAS / COMPAS

Discontinuous radius evolution #978