evaluation based on phase_of_evolution

[mjoyceGR]

Hi all, more phase_of_evolution content:
I ran into serious issues trying to use s% phase_of_evolution >= 2 as a condition for activating a particular modification within other_MLT_results in run_star_extras (this is related to my contextless post about retry: error in do1_dlnT_dm_eqn 3 in bugs-and-problems on Slack). When running an 0.2Ms model at solar composition,

this (case 1) runs fine:

         if (s% phase_of_evolution >= 3) then
            gradr_spot = gradr/factor
         else
            gradr_spot = gradr
         end if

as does this (case 2):

if (s% star_age >= 1d0) then   ...      

but this (case 3) does not:

        if (s% phase_of_evolution >= 2) then
           gradr_spot = gradr/factor
        else
           gradr_spot = gradr
        end if

even though s% star_age = 1 ! yr is definitely after the onset of phase_of_evolution = 2 (the transition time for this parameter combination between phase 1 and phase 2 is about 1d-5 yr)
error.log

"factor" is a number, gradr_spot is the same type as MLT's standard gradr ( type(auto_diff_real_star_order1) )

The result of running case (3) is attached, but can be summarized by
retry: error in do1_dlnT_dm_eqn a bunch of times until failed in do_relax_num_steps.

I confirmed that phase_of_evolution = 1 is met before immediately crashing at phase_of_evolution = 2. I suspect that things work exactly one timestep after the phase_of_evolution = 2 switch comes on, but not at the switch. Tagging @adamjermyn at his request

[rjfarmer]

What's the approximate size of factor?

[rjfarmer]

Seems like is triggering during the create prems relax call. You need to add logic to stop that as the create prems relaxation needs to run with as simple physics as possible to make the initial model.

[mjoyceGR]

The size of 'factor' can vary a lot, but even if I set it to unity, the problem still happens. Hence, I think it is an order-of-evaluations problem rather than a "weird value" problem.

I agree that this appears to be triggered during the relaxation phase, which should be equivalent to phase_of_evolution = 1 (not 2). This issue makes me concerned that phase_of_evolution is not being set at the appropriate time during the transition from relaxation to "real" evolution; namely, that phase_of_evolution switches from 1 to 2 before the mixing routines like other_mlt() have enough information to permit modifying grad_r or grad_a.

Does that make sense?

[rjfarmer]

My point is you shouldn't try to do make any changes to the stellar structure during create_pre_main_sequence. You can check if your in the pre-ms relax phase with (if (s%doing_relax) "relaxing") and I would suggest just returning from other_mlt at that point.

The phase maybe wrong during the relax process, but remember a relax routine is taking the star between point A to point B but there's no guarantee that any intermediate step maps to a real star. During the relax process you should consider it in a "not-a-phase" phase.

[warrickball]

I agree that this appears to be triggered during the relaxation phase, which should be equivalent to phase_of_evolution = 1 (not 2).

mesa/star/private/star_utils.f90

Lines 2794 to 2800 in 680a2b9

	else if (s% L_nuc_burn_total >= s% L_phot*s% Lnuc_div_L_zams_limit) then
	s% phase_of_evolution = phase_ZAMS
	else if (s% log_center_temperature > 5d0) then
	s% phase_of_evolution = phase_PreMS
	else
	s% phase_of_evolution = phase_starting
	end if

I don't think that's what the code does. phase_starting and phase_PreMS are 1 and 2, respectively. I don't see any logic that handles that the star is in a relax phase but that could be added if any relax operations should correspond to phase_starting or a new phase_relax or something.

I agree with @rjfarmer that a reasonable workaround is to handle the relaxation phase in your other subroutine.

PS: Do we ever get models with central temperatures below 10⁵ K...?

[evbauer]

Maybe at the top of this block we should add something like

if(s% doing_relax) then
   s% phase_of_evolution = -1
else if
...

to indicate that we're in a "not-a-phase" state?

[rjfarmer]

Planets might be below 10^5k?

[adamjermyn]

Thanks @mjoyceGR for writing this up.

@evbauer I like that idea. We could create a new phase_relax (which could be -1) and have that get set overriding everything else if we're still in the relax routines.

[evbauer]

Ok, I added phase_relax = -1 in 91047bc.

[mjoyceGR]

Can I request also that phase_of_evolution become a history_columns parameter?

…

On Tue, Nov 9, 2021 at 4:09 PM Evan Bauer ***@***.***> wrote: Ok, I added phase_relax = -1 in 91047bc <91047bc> . — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#340 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADT342O57NJP3R5ZLGTUQDTULGZ4PANCNFSM5HKILPGA> .

-- Dr Meridith Joyce Lasker Fellow Space Telescope Science Institute

[rjfarmer]

I already added it a few days ago in e6055bf

[mjoyceGR]

As discussed with @evbauer last night and with help from @aarondotter and Jamie Tayar (at KITP), I've got a new definition for phase_of_evolution =2 working within run_star_extras. This uses the deuterium birth line (DBL) as the definition for the pre-MS, which is a physically motivated choice unlike the current EEP definition. This is also based on discussions with Marc Pinsonneault, much of which was covered in his 2019 summer school lectures.

See figures

I propose making this the internal definition of phase 2

[warrickball]

Does this work if deuterium isn't in the net?

[mjoyceGR]

probably not, but checks are performed to confirm that deuterium is in the net before trying to determine the physical state. We can just default back to the original phase 2 definition if it's not

[rjfarmer]

What's the new code to determine phase 2? I don't think anyone is really stuck with the old way given its not even been in a release yet, but we should make sure a model without H2 gives a reasonablely similar result (but does not need to be exactly the same).

[mjoyceGR]

still working on it...found some additional issue with Aaron earlier today. I think @evbauer has agreed to help me put this in while at KITP once we're all convinced. I will post plots when I have them.

A few not-minor side issues this also dredges:
(1) having the default setting be pre_ms_relax_to_start_radiative_core = .true. when create_pre_main_sequence_model = .true. is dangerous...this causes a solar-like model to exit the relaxation phase at a central temperature of ~10^6.5, when it should start at something like ~10^5. This will certainly mess up the phase_of_evolution sequence, but I imagine there are many more severe consequences, too. Why was this choice made?

(2) it is strange (to the low-mass star community, from which I hail) that our idea of a "basic" nuclear reaction network does not include deuterium, and likewise that an initial deuterium abundance is not set by default. These are logistical necessities for phase_of_evolution = 2 to work, but if we have a physically motivated definition of the onset of the pre-main sequence, why aren't we using it? I remember talking to @orlox several months ago about what "age zero" actually means. The definition seemed rather arbitrary then (as it does now), but I don't think it has to be.

The reference I'm using is the deuterium birth line (DBL), as described near Figure 1 here: https://ui.adsabs.harvard.edu/abs/2020ApJ...891...29S/abstract

edit: quoted some numbers incorrectly

[rjfarmer]

1. I think the choice was made for massive stars. If we can come up with a criteria for when to switch with low mass stars, then we can add that.

2. Because, from an network standpoint, h2 is only needed to go from h1-> he3 and afterwards its not in any reactions. So we can just make a fake reaction that does h1->h2->he3 and not need to carry the extra isotope. There are always tradeoffs in choosing the network, and basic.net the tradeoff is to have the smallest number of isotopes while getting most of the energy right.

basic.net is really basic and honestly most people should just go to approx21 (but that also lacks h2) so doesn't help here.

If you care about h2 (or any isotope) then you have to pick a network that includes that isotope (and everything that makes/destroys it). If you want we can try to come up with a "low_mass_star.net" network that has just the isotopes that low mass stars care about?

[warrickball]

A few not-minor side issues this also dredges: (1) having the default setting be pre_ms_relax_to_start_radiative_core = .true. when create_pre_main_sequence_model = .true. is dangerous...this causes a solar-like model to exit the relaxation phase at a central temperature of ~10^6.5, when it should start at something like ~10^5. This will certainly mess up the phase_of_evolution sequence, but I imagine there are many more severe consequences, too. Why was this choice made?

From what I remember in discussing the new control with @orlox and Bill, Bill implemented it because many of the changes to things like setting the initial (i.e. near-ZAMS) rotation rate or changing the atmosphere model are much more robustly made after the radiative core has developed in a massive star. That said, I'm in favour of it the default being .false. because it's non-standard and I think it only really matters if you start doing other non-standard things (like including rotation or changing the atmosphere model). Alternatively, we can change some of the default parameters for that part. e.g. I noticed

    pre_ms_check_radiative_core_Lnuc_div_L_limit = 0.1d0

which I think should be much smaller, like 1d-6 or something.

(2) it is strange (to the low-mass star community, from which I hail) that our idea of a "basic" nuclear reaction network does not include deuterium, and likewise that an initial deuterium abundance is not set by default. These are logistical necessities for phase_of_evolution = 2 to work, but if we have a physically motivated definition of the onset of the pre-main sequence, why aren't we using it? I remember talking to @orlox several months ago about what "age zero" actually means. The definition seemed rather arbitrary then (as it does now), but I don't think it has to be.

I think in many circumstances it's reasonable to assume that h2, c13, o15 and other minor intermediate species are rapidly driven to an equilibrium value that can be used instead of tracking that abundance. If you're interested in the details of these abundances, yes, you should change your net.

I'm also not sure what you mean by the "onset of the pre-main-sequence". What comes before the pre-MS? For everything before central hydrogen burning starts is pre-MS.

basic.net is really basic and honestly most people should just go to approx21 (but that also lacks h2) so doesn't help here.

I seldom evolve models that far but don't the defaults automatically switch to approx21 if the evolution goes that far?

[rjfarmer]

which I think should be much smaller, like 1d-6 or something.

Would you mind trying that in a branch and we can see if anything breaks? and if its better for the low mass stars.

I seldom evolve models that far but don't the defaults automatically switch to approx21 if the evolution goes that far?

Yes but net swaps are messy and things are much simpler if you just start with approx21 rather than doing a swap (if you didn't realize you going to do a net swap, you can break your history output when it tries to add new columns for the new isotopes). My point though really is just basic is so basic you should try to avoid it whenever possible.

[warrickball]

which I think should be much smaller, like 1d-6 or something.

Would you mind trying that in a branch and we can see if anything breaks? and if its better for the low mass stars.

See branch whb/pre-ms-Lnuc-limit, where I changed it to 1d-6. Tests look pretty good but six have yet to finish.

[rjfarmer]

@mjoyceGR does @warrickball branch whb/pre-ms-Lnuc-limit give you a better starting temperature for your low mass stars?

[mjoyceGR]

will investigate when there's a break in conference activities

…

On Mon, Nov 15, 2021 at 8:00 AM Robert Farmer ***@***.***> wrote: @mjoyceGR <https://github.com/mjoyceGR> does @warrickball <https://github.com/warrickball> branch whb/pre-ms-Lnuc-limit give you a better starting temperature of your low mass stars? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#340 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADT342MRNH54WVNZP4PSTVLUMEVBPANCNFSM5HKILPGA> .

-- Dr Meridith Joyce Lasker Fellow Space Telescope Science Institute

[mjoyceGR]

(which is probably going to be next week some time...apologies)

[mjoyceGR]

OK, so this could simply be that I'm not working with the branches correctly, but here is what happens

When I run Warrick's branch with

    create_pre_main_sequence_model = .true.
    pre_ms_relax_to_start_radiative_core = .false. 

I get the yellow track, which is fine.

When I run Warrick's branch with

    create_pre_main_sequence_model = .true.
    pre_ms_relax_to_start_radiative_core = .true. 

I get the red (non-)track and this error as soon as the post-relaxation phase starts:

retry: atm get_table: L < 0       1
                           1st model retry log10(dt/yr)   -5.3010299956639804D+00
 retry: max residual jumped -- give up in solver       2
 retry: max residual jumped -- give up in solver       5
 retry: max residual jumped -- give up in solver       5
 retry: max residual jumped -- give up in solver       5

So not only is this starting temperature still inappropriately high, but the atmosphere tables are not happy...Maybe Warrick can share his inlist to make sure I'm not introducing some other kind of problem with my (science) inlist? I'm using photosphere tables, A09, and setting an initial composition. But, all of these things work fine (in the sense of not throwing atm get_table errors) on branch b9a5f96 with create_prems and relax_to_radiative_core both on

[warrickball]

The only thing my branch changes is

pre_ms_check_radiative_core_Lnuc_div_L_limit = 0.1d0

to

pre_ms_check_radiative_core_Lnuc_div_L_limit = 1d-6

in &star_job. What happens if you go to your working branch and change only this parameter in your otherwise working inlist?

[mjoyceGR]

the same thing (which is good, probably):

 retry: atm get_table: L < 0       1
                           1st model retry log10(dt/yr)   -5.3010299956639804D+00
 retry: max residual jumped -- give up in solver       2
 retry: max residual jumped -- give up in solver       5
...
 retry: max residual jumped -- give up in solver       5
 retry: max residual jumped -- give up in solver       5
 retry: max residual jumped -- give up in solver       5
terminated evolution: hydro_failed

[warrickball]

Could you send your inlists (via whatever channel you like) and let me know exactly which commit you're using?

[mjoyceGR]

will do on Slack (don't want to share publicly yet because there are other science components I'm developing)

[warrickball]

Though something of a tangent to the original point of this thread, my cut on pre_ms_check_radiative_core_Lnuc_div_L_limit is probably too restrictive at 1d-6. But I've spun this off into a separate issue (#347) because it's tangential to this one.

[mjoyceGR]

FYI, I have made a branch on which I am incorporating the DBL definition for phase_of_evolution = 2with @evbauer's help. Still figuring out how to get this on github, but it's coming.

Something else worth noting here and in #347 is that pre_ms_T_c is set to a default value of 3d5 regardless of whether pre_ms_relax_to_start_radiative_core is being used, but both the DBL and non-DBL definitions of phase_PreMS tests against log_center_temperature > 5. So, every star is being "born" with a higher core temperature (by default) than this...I think that is a bad default configuration.