Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated wind perscriptions #996

Merged
merged 42 commits into from
Nov 2, 2023
Merged

Updated wind perscriptions #996

merged 42 commits into from
Nov 2, 2023

Conversation

jmerritt1
Copy link
Collaborator

@jmerritt1 jmerritt1 commented Sep 27, 2023

Added a new UPDATED option to the mass-loss prescriptions. This is a wrapper applying new specific prescriptions on a case-by-case basis. These are selected via new program options for OB, VMS, RSG, and WR-mass-loss, which each have their own wrapper function selecting from a total of 15 new perscriptions. Previous default behavior is recoverable by specifying VINK.

Commit history is in the updated_winds branch, which became relatively unmergable. The changes have been collected and added to an up-to-date branch here as Jeff suggested.

…er merge, changes broken down with history in updated_winds
@SimonStevenson
Copy link
Collaborator

SimonStevenson commented Sep 27, 2023

Two things that will need updating:

  • Update changelog.h with description of changes, and bump version
  • Update documentation for new options

@SimonStevenson SimonStevenson added the enhancement New feature or request label Sep 27, 2023
@jmerritt1
Copy link
Collaborator Author

note in Options.cpp, that we didn't use, (" + AllowedOptionValuesFormatted(" and instead list the options as before. This convention must have been changed in the last few months. See

("OB mass loss prescription (options: [NONE, VINK2001, VINK2021, BJORKLUND2022, KRTICKA2018], default = " + p_Options->m_OBMassLoss.typeString + ")").c_str()
for example.

@jeffriley
Copy link
Collaborator

note in Options.cpp, that we didn't use, (" + AllowedOptionValuesFormatted(" and instead list the options as before. This convention must have been changed in the last few months. See

("OB mass loss prescription (options: [NONE, VINK2001, VINK2021, BJORKLUND2022, KRTICKA2018], default = " + p_Options->m_OBMassLoss.typeString + ")").c_str()

for example.

Yeah, that changed a while ago - we should use the new method in this PR (not difficult)

@jeffriley
Copy link
Collaborator

missed one new RSG option on ln 1736, but should change to AllowedOptionValuesFormatted() anyway?

Yes - that way the help text is automatically updated if any of the allowed values are changed - and the help text format is kept uniform with other options.

Jonathan Merritt added 2 commits September 27, 2023 16:44
@jmerritt1
Copy link
Collaborator Author

missed one new RSG option on ln 1736, but should change to AllowedOptionValuesFormatted() anyway?

Yes - that way the help text is automatically updated if any of the allowed values are changed - and the help text format is kept uniform with other options.

Thanks, one less spot to change options in the future. I've updated these.

@ilyamandel
Copy link
Collaborator

Sorry, I am offline until October 8, so won't be able to review this -- defer to @jeffriley and @SimonStevenson !

@jmerritt1
Copy link
Collaborator Author

jmerritt1 commented Sep 28, 2023

Two things that will need updating:

  • Update changelog.h with description of changes, and bump version
  • Update documentation for new options

in program-options-list-defaults.rst, it's formatting some of the options lists with the bracket on the next line, but I'm not seeing why?

@jeffriley
Copy link
Collaborator

bumped version to 02.39.02

Let's make it 02.40.00 - this is new functionality, not a defect repair

xx.yy.zz -> xx is (supposed to be) major new release; yy minor new release/new functionality; zz defect repair

@jmerritt1
Copy link
Collaborator Author

bumped version to 02.39.02

Let's make it 02.40.00 - this is new functionality, not a defect repair

xx.yy.zz -> xx is (supposed to be) major new release; yy minor new release/new functionality; zz defect repair

Should I list all the specific perscriptions in the changelog or are just the new program options fine?

Copy link
Collaborator

@jeffriley jeffriley left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @jmerritt1. I have a few comments and suggested changes. I'll number them so we can refer to them by number if necessary:

  1. prescription is mis-typed as perscription in a few places:

    program-options-list-defaults.rst
    BaseStar.cpp
    changelog.h
    constants.h
    Options.h

  2. Should we implement an option to set the VERY_MASSIVE_MINIMUM_MASS threshold, or is 100Msol either the universally accepted threshold or required for our implementation?

  3. If you run your version of the executable with ./compas --help, you will see:

    Options:
    ...
    --OB-mass-loss
    OB mass loss prescription (Options: [], default = 'NONE)
    ...

    Note no options listed in [] - that's because you missed step 9 in the instructions at the top of Options.cpp

    Also note that there is a missing single quote after 'NONE - that needs to be added in Options::AddOptions():

    Change

         `("OB mass loss prescription (" + AllowedOptionValuesFormatted("OB-mass-loss") + ", default = '" + p_Options->m_OBMassLoss.typeString + ")").c_str()`
    

    to

        `("OB mass loss prescription (" + AllowedOptionValuesFormatted("OB-mass-loss") + ", default = '" + p_Options->m_OBMassLoss.typeString + "')").c_str()`
    

    and similarly for RSG-mass-loss, VMS-mass-loss, and WR-mass-loss

  4. If you run your version of the executable with ./compas --create-yaml-file myyaml, and look at the resultant yaml file (myyaml). At the bottom of the file you'll see:

     ### Additional COMPAS options not found in YAML template ###
     #    --RSG-mass-loss: 'NJ90'                                               # Default: 'NJ90'
     #    --VMS-mass-loss: 'NONE'                                               # Default: 'NONE'
     #    --WR-mass-loss: 'BELCZYNSKI2010'                                      # Default: 'BELCZYNSKI2010'
     #    --metallicity-distribution: 'ZSOLAR'                                  # Default: 'ZSOLAR'          Options: ['LOGUNIFORM','ZSOLAR']
    

    That's because your changes to yaml.h:

         "    --OB-mass-loss: 'NONE'                                                # Default: 'NONE'"
         "    --RSG-mass-loss: 'NJ90'                                               # Default: 'NJ90'"
         "    --VMS-mass-loss: 'NONE'                                               # Default: 'NONE'"
         "    --WR-mass-loss: 'BELCZYNSKI2010'                                      # Default: 'BELCZYNSKI2010'"
    

    are malformed. You need to remove the default values in both places. e.g. change

         "    --OB-mass-loss: 'NONE'                                                # Default: 'NONE'"
    

    to

         "    --OB-mass-loss"
    

    and similarly for RSG-mass-loss, VMS-mass-loss, and WR-mass-loss (see examples above and below your additions in yaml.h).

  5. BaseStar::CalculateMassLossRateOBBjorklund2022():

    • function name in description does not match actual function name
    • is it worth noting the caveats documented in the description (i.e. "This prescription is calibrated to the following ranges:") in the online documentation for the relevant option? Is it important a user knows that (particularly a user who may just use COMPAS without looking at the code)?
  6. BaseStar::CalculateMassLossRateOBVinkSander2021():

    • function name in description does not match actual function name
    • function has commented code that should be removed
    • equations at lines 2015-2020, 2030-2035, and 2045-2051 have ragged right edges with math operators. Either have a single space between the end of the text and the operator, or line all operators up (like the equation at 1967-1973) - the way it is makes it too hard to read.
    • it's easier to read if the comments at the end of code lines are aligned with each other too (for humans maintaining the code - the compiler doesn't care of course...)
  7. BaseStar::CalculateMassLossRateOBKrticka2018():

    • function name in description does not match actual function name
    • log10(m_Metallicity / ZSOL) is BaseStar::LogMetallicityXi() - log10() is expensive, no need to call it if we don't have to (and twice...).
  8. BaseStar::CalculateMassLossRateRSGBeasor2020() has commented code that should be removed.

  9. BaseStar::CalculateMassLossRateRSGDecin2023():

    • I would probably just rewrite this as:

       return PPOW(10.0, -20.63 - 0.16 * m_MZAMS + 3.47 * log10(m_Luminosity));
      

      but, can't we rewrite that as:

       return PPOW(m_Luminosity, -20.63 - 0.16 * m_MZAMS + 3.47);
      

      I think the latter would use fewer compute cycles - the former seems to do an unnecessary log10() operation, which isn't cheap.

  10. BaseStar::CalculateMassLossRateRSGKee2021():

    • has commented code that should be removed.
    • also, change 2. -> 2.0, 17000. -> 17000.0, and 60000. -> 60000.0 (for the same reason as above - easier for humans to read (esp. at a glance))
  11. BaseStar::CalculateMassLossRateRSGVinkSabhahit2023():

    • change -8. -> -8.0, and -24. -> -24.0 (same as above)
  12. BaseStar::CalculateMassLossRateVMSBestenlehner2020():

    • change 1. -> 1.0 (same as above)
    • change constants 1 and 2 to 1.0 and 2.0. c++ is a strongly-typed language: "1" is an integer, whereas "1.0" is a double. It may make no difference now, but it might if someone changes things later - it's good practice to use the correct data type even for constants.
  13. BaseStar::CalculateMassLossRateVMSVink2011():

    • reference/link to paper in the description (is 'Vink 2011' sufficient)?
    • function should be const
    • if the else clause is executed, rate is calculated, and is recalculated after the if/else using logMdotdiff which is undefined (only defined in the if clause is executed). teff and rate2001 only need to be calculated in the else clause, so those calculations should be moved into the else clause. Bring the final calculation of rate up into the if clause (it requires logMdotdiff).
  14. BaseStar::CalculateMassLossRateVMSSabhahit2023():

    • reference/link to paper in the description (is 'Sabhahit 2023' sufficient)?
    • function should be const
    • resolve and remove question on line 2260
  15. BaseStar::CalculateMassLossRateWolfRayetSanderVink2020():

    • change the conditional on line 2425 to use utils::Compare()
    • Mdot has already been initialised to 0.0 (at line 2412), so no need to set it 0.0 at line 2426. Instead, change the conditional to
      if (utils::Compare(logL0, logL) < 0) and delete lines 2426 and 2427.
    • move line 2434 into the if - no need to multiply by the option value if Mdot is 0.0 - just wastes compute cycles.
    • change - 1 to - 1.0 on line 2430
  16. BaseStar::CalculateMassLossRateVMS() has commented code that should be removed

  17. BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023():

    • change 141E3 to 141.0E3 on line 2457 (same as above re data type)
    • change 100E3 to 100.0E3 on line 2458 (same as above re data type)
    • change the conditional on line 2465 to use Utils::Compare()
  18. BaseStar::CalculateMassLossRateHeliumStarVink2017(): I think I'd just use a single return line, but up to you.

  19. BaseStar::CalculateMassLossRateHeliumStarShenar2019() has commented code that should be removed

  20. BaseStar::CalculateMassLossRateOB(): parameter should be const

  21. BaseStar::CalculateMassLossRateRSG(): parameter should be const; function should be const

  22. BaseStar::CalculateMassLossRateVMS(): function should be const

  23. BaseStar::CalculateMassLossRateUpdatedPrescription():

    • "The structure is similar to above" - do you mean similar to BaseStar::CalculateMassLossRateVink()? If so, say so (it may not be above if the code gets changed/refactored).
  24. Add an entry to the What's New page of the online documentation to describe any changes to default behaviour (e.g. which mass loss prescription is now default? Any other changes to default behavour?)

  25. Your question: "Should I list all the specific perscriptions in the changelog or are just the new program options fine?"
    - Just the new options are fine for changelog, but you should add a more complete descripion in the 'What's New' entry (see above).

  26. Given the number of times we use log10(m_Luminsoity), I wonder if it's worth calculating that when m_Luminosity is calculated, and storing it in a class member variable (call it m_Log10Luminosity - similar to m_Log10Metallicity)?

That's about it - there are some minor formatting changes I'll either point out later or just fix myself at the final review.

I have not checked the functionality - I checked that the code runs and data is produced, but I have not checked if the contents of the data files are correct according to the new prescriptions added, or whether any existing functionality has been negatively impacted (wrt the correctness of the data produced). I assume you have satisfied yourself that the data produced is correct.

@jeffriley
Copy link
Collaborator

And...

in program-options-list-defaults.rst, it's formatting some of the options lists with the bracket on the next line, but I'm not seeing why?

I think that's just the way Read the Docs formats - the line break just happens to be there. You could play around with the text to shift the line break away from the bracket if it really bothers you.

@jeffriley
Copy link
Collaborator

One more change:

Change

constexpr double LSOLW = 4E26;

in constants.h to

constexpr double LSOLW = 4.0E26;

@jeffriley
Copy link
Collaborator

Wait...

but, can't we rewrite that as:

return PPOW(m_Luminosity, -20.63 - 0.16 * m_MZAMS + 3.47);

No, that's wrong (I mentally put parentheses around -20.63 - 0.16 * m_MZAMS + 3.47). Can we still avoid the log10? Maybe...

…g file (dont use the method paper config)"

This reverts commit 346e98d.
@SimonStevenson
Copy link
Collaborator

SimonStevenson commented Oct 31, 2023

Thanks everyone for the input on this PR, it's great that it's almost there.

running with the new default prescription won't artificially reduce WR by an extra factor of 10?

Currently, it does reduce WR mass-loss rates by an additional factor of 10. I believe that all of JDs testing/results have been performed with this setup. I noticed this the other day. We will do some additional testing/comparisons of fWR=0.1 v fWR=1.0. I think that the new WR mass-loss prescription does not reduce mass-loss rates by as much as a factor of 10 for most high-mass He stars.

I made a plot comparing the initial-final He star masses (at solar metallicity) for both the old and new prescriptions, with fWR=0.1 and fWR=1.0.
heliumStarInitialFinalMassLossRelation

Here we see that whilst for low-mass He stars (where we use Vink2017) the mass-loss rates are reduced, for high-mass solar metallicity He stars we actually see increased mass loss.

running with the BELCZYNSKI2010 prescription will reproduce the immediately preceding default (including that factor of 0.1)?

Currently this is not the case -- whatever value of fWR is set will be used. If the user doesn't specify a value, the default will be used (currently 0.1).

Based on this, I'm not sure what the new default fWR should be (fWR=1.0 better reflects the implemented mass-loss prescriptions, fWR=0.1 allows for the production of Cyg X-1 like BH masses at solar metallicity).

@jmerritt1
Copy link
Collaborator Author

jmerritt1 commented Oct 31, 2023

Thanks everyone for your work on this. I've been travelling for two weeks but just got back online and caught up with the discussion.

A. @jmerritt1, from your comment above:

... m_DominantMassLossRate ... It is my preference to only set this flag in higher wrapper functions ...
Can you just do a final check that m_DominantMassLossRate gets set when it should (e.g. scenario above - BaseStar::CalculateMassLossRateOBVink2001() is called from multiple places).

I've checked again that m_DominantMassLossRate is set in both the BELCZYNSKI2010 (added one missing line for OB) and FLEXIBLE2023 wrappers. For some of the giant branch specific prescriptions, it is set in the relevant stellar type files under the Hurley wrapper.

B. enum class MASS_LOSS_TYPE in constants.h somewhat confusingly mixes mass loss as applied to different types of stars (VERY_MASSIVE, RED_SUPERGIANT, LBV) with specific prescriptions (say, VINK -- Vink wrote a lot of papers about mass loss from different kinds of stars, so not clear what this refers to in this context)

I changed this list to include broader types of mass loss, and not specific prescriptions. We are leaning towards only including types, eg. RSG, VMS, etc. I agree that VINK should be removed, since it should now no longer be set in the higher level wrappers in favor of OB. NIEUWENHUIJZEN_DE_JAGER, KUDRITZKI_REIMERS, and VASSILIADIS_WOOD are sometimes set by the stellar type specific giant branch files. If distinguishing between these isn't an important diagnostic tool, we can change all these to RSG. Is there any reason we need the older prescription specific labels?

Regarding C. and fWR, I've been running with fWR=1.0 for everything so far, and hadn't realized the default .yaml had changed.

@ilyamandel
Copy link
Collaborator

@SimonStevenson -- I am trying to match up your latest plot with the one from your paper with Teagan. That one was in units of luminosity rather than mass on the abscissa -- but it seemed to show that at the highest luminosities (i.e., high masses), mass loss rates from newer models were actually below those from the old models with a rescaling factor of 0.1, right?

image

@ilyamandel
Copy link
Collaborator

@jmerritt1 -- I like your suggestion of removing VINK and replacing all of NIEUWENHUIJZEN_DE_JAGER, KUDRITZKI_REIMERS, and VASSILIADIS_WOOD with RSG. Thanks for doing that!

@SimonStevenson
Copy link
Collaborator

SimonStevenson commented Oct 31, 2023

@ilyamandel We're not using Vink2017 for high-mass (high luminosity) WR stars, only for low-mass (low-luminosity) stripped He stars. At solar metallicity, the transition is around L = 10^5 Lsol (which corresponds to an initial helium star mass somewhere in the range 5-10 Msol). We are using the prescription from Sander & Vink 2020 for high-mass He stars. This results in similar mass-loss rates to the old Belczynski+2010 prescription.
compareVink2017SanderVink2020Solar

@ilyamandel
Copy link
Collaborator

Thanks for the explanation, @SimonStevenson ! I wonder if we ultimately want two separate WR multiplier choices, for Vink 2017 at lower masses and for Sander & Vink 2020 at higher masses? If we did, I might be tempted to use fWRlowmass = 1 and fWRhighmass=0.1 in my runs... But I think the default should be 1 for both, since these are the best available models, and so long as we use just one multiplier, let's keep the default at 1.

BTW, does the transition happen exactly at the point where the Vink 2017 and Sander & Vink 2020 curves cross over in your plot? If not, wind mass loss rates would be discontinuous, which is a bit weird.

@SimonStevenson
Copy link
Collaborator

Yes, we just use the maximum of the two prescriptions, so the transition is where they cross over.

@ilyamandel
Copy link
Collaborator

OK, then my only two remaining requests are:

  • remove VINK and replace all of NIEUWENHUIJZEN_DE_JAGER, KUDRITZKI_REIMERS, and VASSILIADIS_WOOD with RSG in MASS_LOSS_TYPE (and in the code where it's assigned)
  • set the default Wolf-Rayet multiplier to 1.0 rather than 0.1 (and regenerate the default .yaml)

@jmerritt1
Copy link
Collaborator Author

jmerritt1 commented Oct 31, 2023

@jmerritt1 -- I like your suggestion of removing VINK and replacing all of NIEUWENHUIJZEN_DE_JAGER, KUDRITZKI_REIMERS, and VASSILIADIS_WOOD with RSG. Thanks for doing that!

Ah actually I've forgotten we will need to add a label for evolved/giant stars that are also below the cutoffs for RSG, see the pink lines (KR) in the below figure for example. (colored by dominant mass loss rate)
Also, is there any opposition to using acronyms, RSG, LBV, in place of full labels, since this is mainly a diagnostic tool?

image

@jmerritt1
Copy link
Collaborator Author

jmerritt1 commented Oct 31, 2023

And in the above figure, just want to make sure that this is the intended treatment for the gap between 12.5kK and 8kK, around HG, Simon and I came to the understanding that you'd like these to be treated as WR (see short purple line segments around HG)? We revisited the original conversation with Andreas Sander and this seems to be the intention?

@ilyamandel
Copy link
Collaborator

ilyamandel commented Nov 1, 2023 via email

@jeffriley
Copy link
Collaborator

and regenerate the default .yaml

Just in case you missed it, there's now an option to do that (--create-yaml-file, see https://compas.readthedocs.io/en/latest/pages/User%20guide/Program%20options/program-options-list-defaults.html#options-props-c)

@jmerritt1
Copy link
Collaborator Author

jmerritt1 commented Nov 1, 2023

and regenerate the default .yaml

Just in case you missed it, there's now an option to do that (--create-yaml-file, see https://compas.readthedocs.io/en/latest/pages/User%20guide/Program%20options/program-options-list-defaults.html#options-props-c)

Just tried this but have it writing the previous default, fWR=0.1. Where should I change this? options.cpp? (edit, just wasn't working because I forgot to recompile before generating).

@jeffriley
Copy link
Collaborator

jeffriley commented Nov 1, 2023

@jmerritt1 I don't think you pushed the new version of Options.cpp - with m_WolfRayetFactor = 1.0

Oh - and compasDefaultConfig.yaml should be in compas_python_utils/preprocessing, not src

@jmerritt1
Copy link
Collaborator Author

@jmerritt1 I don't think you pushed the new version of Options.cpp - with m_WolfRayetFactor = 1.0

Oh - and compasDefaultConfig.yaml should be in compas_python_utils/preprocessing, not src

Okay I think it's fixed now, and I've removed the one in src.

Copy link
Collaborator

@ilyamandel ilyamandel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good -- thanks for all the hard work, everyone!

@ilyamandel ilyamandel merged commit e39c6b7 into dev Nov 2, 2023
2 checks passed
@ilyamandel ilyamandel deleted the updated_winds_fresh branch November 2, 2023 07:44
@ilyamandel
Copy link
Collaborator

@jmerritt1 -- I saw your new issue post via e-mail, but it seems to have disappeared. I suggest posting new issues separately, since people will likely miss them in closed PRs.

jeffriley pushed a commit to jeffriley/COMPAS-1 that referenced this pull request Dec 15, 2023
reinhold-willcox added a commit that referenced this pull request Dec 22, 2023
* Fix for issue #1022

* Code cleanup

* Event strings (#1021)

* (compasUtils) bug fix to getSNevents()

* (compasUtils) make getEventHistory() work for COMPAS output with unordered seeds

* (compasUtils) fix typo

* switched to lexsort from argsort to better handle the ordering of multiple events from the same seed

* changed list initialization for np array, changed some variable names for clarity

* changed np array back to list, since array couldn't handle input of different sizes

---------

Co-authored-by: Reinhold Willcox <[email protected]>

* Merge pull request #996 from TeamCOMPAS/updated_winds_fresh

Updated wind perscriptions

* Reapply winds changes from 2.40.00

* fixes to the detailed plotter, including improved robustness to the Roche Lobe, and masking for only end-of-timestep events

* added question to CalculateMassLossRateNieuwenhuijzenDeJager

---------

Co-authored-by: Jeff Riley <[email protected]>
Co-authored-by: Ilya Mandel <[email protected]>
Co-authored-by: Mike Lau <[email protected]>
Co-authored-by: jeffriley <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants