Added new rerun type: dedt_energy_eqn

[sgossage]

This rerun is designed to help with certain numerical issues experienced by stars undergoing superthermal accretion and pulsations during late stage stripped He star evolution.

This rerun (dedt_energy_eqn) will...

1. Enable MESA's dedt form of its energy equation. This helps numerical convergence during rapid (superthermal) accretion.

Three new termination codes are introduced as well to catch cases of contact binary evolution when the accretor is inflated due to superthermal accretion (one termination code) and when the accretor is rotating at a critical rate while undergoing superthermal accretion (two termination codes for the depending on which star is the current donor). The former is a contact binary case, and the latter are akin to L2 overflow cases (loosely based on physical arguments from e.g., Lu+ 2023).

Most of the changes take place in run_binary_extra.f90. The associated changes may be found here: POSYDON-MESA-INLISTS/commits/seth_dedt_energy_eqn/

This also fixes an issue where the current wind prescription was not stored as a unique value for each star.

Additional changes that may enter after further testing for stripped He stars (WIP)

This primarily helps lower mass ratios where rapid (superthermal) accretion happens. However, stripped He stars occur at higher mass ratios on our HMS-HMS grids, and on the CO-HeMS grid. Those need help beyond changing the energy equation. These may need alterations to MLT++ which are currently being tested. 3. ~~Enable MESA's `lnPgas_flag` when a star forms a helium core. This causes MESA to use Pgas, rather than density to solve the energy equation, which helps numerical stability with stripped He stars.~~ 4. ~~Enable MESA's `v_flag` when a star becomes a stripped He star. This helps numerical convergence during pulsations that occur during late stage stripped He star evolution. Also set `convective_bdy_weight=0` at this point because this algorithm often causes segmentation faults after models start to settle as white dwarfs.~~ 5. ~~Introduces logic to use a boosted (Bloecker at the moment) wind for stripped He stars.~~

Details from previous version of this PR (deprecated)

These rerun types address issues certain models have with numerical convergence; these are described below. - [x] `rerun_type = timestep_MTlimit` (HMS-HMS fix, source code[†](https://github.com/POSYDON-code/POSYDON-MESA-INLISTS/tree/seth_HMS-HMS_rerun_accretion_fix)) - [x] `rerun_type = egrav_Pgas_dXdt_terms` (CO-HeMS fix, inlist tweaks[†](https://github.com/POSYDON-code/POSYDON-MESA-INLISTS/tree/seth_CO-HeMS_rerun_fix)) 1. **`rerun_type = timestep_MTlimit`:** On the HMS-HMS, this tends to occur with an accretor, near L2 overflow conditions. Relatively high mass transfer rates, super adiabatic conditions (w/ He ionization opacity bumps at the surface of both stars), and super-critical rotation of the accretor tend to conspire to cause MESA's hydrodynamics solver to fail to converge on a solution. This limits the accretion rate once the timestep falls below 10 years, and completely turns it off once the timestep falls below 1 day. This is meant to let the accretor star's structure settle and for the timestep to recover. This functionality will be enabled if `x_logical_ctrl(7) = .true.` in either inlist1 or inlist2. [†link to branch](https://github.com/POSYDON-code/POSYDON-MESA-INLISTS/tree/seth_HMS-HMS_rerun_accretion_fix) 6. **`rerun_type = egrav_Pgas_dXdt_terms`:** Perhaps the most numerous of the CO-HeMS models failing to converge experience numerical issues associated with high nuclear burning and superadiabaticity. Switching the basic variable MESA uses to solve `eps_grav` from density to gas pressure aids numerical stability, as does enabling usage of the composition terms in calculating `eps_grav`. The changes in the inlists are ``` include_composition_in_eps_grav = .true. change_initial_lnPgas_flag = .true. new_lnPgas_flag = .true. ``` [†link to branch](https://github.com/POSYDON-code/POSYDON-MESA-INLISTS/tree/seth_CO-HeMS_rerun_fix) Additionally, a number of models that originally experienced segmentation faults are able to converge when `convective_bdy_weight = 0` is used; apparently this algorithm has something like memory alloc. issues in some cases, so disabling it helps. The pre-existing rerun type `conv_bdy_weight` should be used for these. A third class of very low mass (<= 0.7 Msol), very short period (< 0.1 days) systems still fail to converge in the CO-HeMS grid (redbacks?). These tend to evolve further with `use_dedt_form_of_energy_equation = .true.`, but this option affects the evolution of the other models, and even when enabled, not all of these systems are solved. So, this last option is not used for now and these models remain unaccounted for.

[sgossage]

Right now rerun_type = co-hems_MLTpOFF only checks for the 'reach cluster timelimit' termination flag. Is that what most of the CO-HeMS models terminate with? Should this rerun_type check for more flags?

[mkruckow]

About 1., why do you only check for the time step of star 1? Wouldn't it be better to check the minimum for star 1's and star 2's time step is in the window between 1 day and 1 yr?

About 2., I wouldn't include the 'convective_bdy_weight' : 0, where it isn't needed. As you said we'd have a different rerun for this. Another question, I'd have here: the models, which have issues with MLT++, don't they may stop with min_timestep_limit? Hence, is there a good reason to not sure the TF1_POOL_ERROR?

In general: While you named the reruns to be meant for a specific grid type, you didn't check for this while exporting the rerun. Hence, either remove the grid specification in the rerun name, if it could be run on any grid or filter for the grid type when exporting reruns.

[sgossage]

Thanks, Matthias.

1. I've changed this to check for the min timestep between star 1 and 2 now, as suggested. Originally I was going off the thinking that the code already checks for the minimum in $MESA/binary/private/binary_timestep.f90, but may be better to just do it explicitly to not worry about function call order, etc.
2. I've removed 'convective_bdy_weight' : 0 now. From the test models I have, all of the runs fixed by this rerun type terminated originally due to wall time, but I've changed the termination_flag to TF1_POOL_ERROR for now, just in case. Some of the min_timestep cases are fixed typically by setting convective_bdy_weight = 0 and some remain unresolved (those that remain unresolved are helped by use_dedt_form_of_energy_eqn = .true., but it's not a solution for 100% of them).
3. I've removed the grid specification from the rerun names as suggested because I think these could be run in principle on any grid that might need them.

[mkruckow]

1. Yeah, MESA will take the minimum of all the time steps finally to determine the new time step. In the binary module there are three sources for the minimum time step:

   1. star 1
   2. star 2
   3. binary interaction

   Each of the three has a bunch of different checks how to determine the minimum time step. You are changing a parameter of the binary interaction in case the star(s) can't handle it. I think that's fine, but as you already changed it should allow to take care of both stars having issues. Btw. we do have another correction on the parameter b% mass_transfer_beta later in the same file (about line 1500 of run_binary_extras.f). I think it would be better to put those things together and this about, which should overrule the other. One more comment to your recent change in your MESA-inlists branch, you need to put the check whether a star is a point mass back in place, because b% s1 won't be well defined in case s1 is a point mass (the same holds for s2). Hence:

   min_dt = t_hi
   if (b% point_mass_i /= 1) then
   min_dt = min(min_dt, b% s1% dt)
   if (b% point_mass_i /= 2) then
   min_dt = min(min_dt, b% s2% dt)

2. OK, thanks.
3. OK, thanks.

[sgossage]

Thanks again, I've implemented these suggestions.

Could you clarify a bit for your point on code near line 1500? Are you suggesting we would replace that code or retain it and introduce a switch between this PR's functionality and that part of the code?

[mkruckow]

Sure, I can be a bit more detailed: I had two things in mind:

1. move the new code there, to have the modifications on the same variable closer together (I'm not sure whether you have been aware of this code part, when you did your implementation.)
2. Because the code parts do modify the same parameter we need to think about, what should happen if both cases are fulfilling at the same time. Usually, the later code will simply overwrite what was modified before. The already existing part was implemented to keep the mass transfer solution more stable close to critical rotation, hence can you tell, whether this is used in the cases you looked into and why it doesn't help or why the stars in your case don't get to critical rotation before running into other issues. There is one more point, the old code turns the mass_transfer_beta back to 0, in case it is 1 and there is a low mass transfer rate. Hence, are there cases, where the old code would may prevent your code from being effective?

[mkruckow]

Please don't forget to add the new rerun types to the documentation [docs/_source/components-overview/post_processing/pipeline/pipeline_steps.rst].

[sgossage]

Sure, I can be a bit more detailed: I had two things in mind:

1. move the new code there, to have the modifications on the same variable closer together (I'm not sure whether you have been aware of this code part, when you did your implementation.)
2. can you tell, whether this is used in the cases you looked into and why it doesn't help or why the stars in your case don't get to critical rotation before running into other issues. There is one more point, the old code turns the mass_transfer_beta back to 0, in case it is 1 and there is a low mass transfer rate. Hence, are there cases, where the old code would may prevent your code from being effective?

I see, thanks.

RE 1: OK, will do. I think maybe using an x_logical_ctrl could be good to enable/disable this feature. Maybe there are cases it shouldn't be used

RE 2: This isn't used in the cases I've tested with; for these cases x_logical_ctrl(4) = .false.. My impression is that it's false for the rest of the grid too, but is that the case? I agree that the code mentioned could e.g., turn MT back on when the code I'm implementing wants it to be off (it would be turned on if MT rate <= 1e-7 Msol/yr if mass_transfer_beta = 1). That might not be a problem since MT is low, but I think making the new code being introduced conditional on an x_logical_ctrl as well could be a good compromise for now, since utilizing both pieces of code is untested at the moment and unintended (but not necessarily conflicting).

[mkruckow]

I talked with @ZepeiX, who implemented the other part of the code. He confirmed, my check, that we never used his code there, because it was meant to help before the index bug for reverse MT was found in MESA. Thus, the other code is probably out of date and we don't need to worry too much about. But having a control flag to switch your change on/off might be useful.

[mkruckow]

Just another thought I just had: For the reduction of the accretion efficiency, should we request a minimum mass transfer rate here, e.g. 1e-5 Msun/yr? I think this would already go in the direction to more apply it to systems, which do have significant mass transfer hence although having the accretor to have some material being accreted beforehand and therefore potentially being faster rotating as well? Or is this fix needed for slowing rotating stars already on the very early stage of mass transfer?

[mkruckow]

Update (March 28): The second rerun type is ready, there are some tries on using a different function for the first one, hopefully completing till next week.

[astroJeff]

Update (Apr 18): PR consists of links to the git hashes for two new re-runs. The re-run for MT thermal timescale limit is still being perfected. The other re-run (Pgas equations) and not clear how helpful this re-run actually is. It helps at 2 Zsun for HMS+CO RLO grid. Maybe not so much for low Z grids. We will discuss during Monday's collaboration meeting and come up with plan to resolve this.

[astroJeff]

Update (Apr 25): @sgossage is currently testing these re-runs. Will report back next week

[astroJeff]

Update (May 9): @sgossage to update the rerun specifics based on presentation at last collaboration meeting.