Spurious behaviour of retry during mass transfer

[Daniel-Pauli]

Hi everyone,

I recently wanted to update my MESA version to the newest version. When playing around and testing my models to see if they converge properly, I stumbled over some issues with the retry happening during mass-transfer, when the secondary rotates close to critical rotation (typically encountered in my Case B models). Since I am using custom ZAMS etc. I tried to make a simple testcase with the "standard" binary file such that it can be reproduced easily.

The thing is that to get a better resolution during the mass-transfer I used in previous versions the delta_HR_hard_limit (and some other limits, but this is already enough to showcase it) to avoid some weird jumps in the HRD. If I use this in the new version the code also tries to make a retry but then "jumps away" from its original position and leads to smaller and smaller timesteps, until the code crashes. However, if I turn off the hard limit (in the attached files I do a trick and turn if off only during mass transfer (see run_binary_extras), to still be able to resolve the TAMS) the model runs though the mass transfer and only has some small jumps.

I attached a PDF showing an HRD comparing the two models. One can see that the blue track diverges a bit. I have set the min_timestep_limit to ~3, if one would make it smaller it would diverge even more. However, the red dashed line continues and shows some spikes later. I have already checked if the delta_HR_hard_limit is triggered correctly or if it triggers a retry at the wrong place. But if my math is correct, the retry is triggered at the correct position, but leads to strange results.

I am not sure if this is a bug, or linked to some inlist command that I am not aware of. I would be happy to get some feedback.

For those who want to test this, I uploaded my work directory of this "simple" testcase to this post. If you do not want to run the models by yourself, I also uploaded them to my google drive (https://drive.google.com/drive/folders/1S5Cnrghax2O6Ym0fmCyQGDIjBiRDbDPS?usp=drive_link and https://drive.google.com/drive/folders/1_LT21HuyBgBkCCiYXfrFqElN19lWFQZF?usp=drive_link) but I will delete the files in the next 2 weeks. Note, that I saved the terminal output in "out.txt" in the folders, in case you want to have a look on whats going on in the models.

I hope that most of the things I described are clear, if not please let me know and I will elaborate further on this issue.

All the best, Daniel Pauli

HRD.pdf test_retry_during_mtrans.zip test_retry_during_mtrans_no_hard_limit.zip

[Debraheem]

Hello, I tried to run your models, but they fail converge on the first timestep? Are you sure these are the exact inlists you used to produce your "out.txt" file? Also, I did a grep on both folders, and a diff on the inlists, and both model files you shared are identical and neither have an explicit control for "delta_HR_limit" listed in them?

• I used version r23.05.1

Nonetheless, I would reccomend relaxing "delta_HR_limit" if your hard limits are causing the model to crash. You could be inadvertantly preventing the numerical solver from accepting a certain solution that requires the model to move in the HR diagram, and or an initial timestep dive could also be causing you to resolve a dynamical feature like a pulsation, which will crash the model if hydro (v_flag) is off.

Looking at your out.txt file, most of the solver residuals are coming out of the the luminosity equation "max resid equL", which is a known issue that's being worked on and documented in the current MESA version, since the current version of MESA is more stringent with the solver (where previous version were much more relaxed "sloppier" with numerical residuals). See https://github.com/MESAHub/mesa/issues/563.

You could try "convergence_ignore_equL_residuals = .true.", and see what happens.

If you can re-share a working version of your inlists I can try giving it another shot.

[Daniel-Pauli]

Hi Eb,

thanks for your help. I am also using r23.05.1. I just noticed that I might be using a custom burning net... Sorry about that! If you want to change it, just use "basic.net" and "co_burn.net" and then it should work. I will try your suggestion of more relaxed solver residuals. I will come back to this when my models are finished!

About the inlists, indeed the delta_HR_limits are set the same (also the hard limit) I did not noticed that. But I am turning delta_HR_hard_limit on and off in the run_binary_extras.f90 after mass-transfer started, because I wanted to have still the high resolution at the TAMS and that I am comparing apples with apples, as I force the model to behave similarly until the onset of mass transfer. For more details see extras_binary_finish_step contains

         !remove delta_HR_limit_hard during mass transfer because a retry leads to spurious behavior
         if(abs(b% mtransfer_rate) > 1d-10) then
            b% s1% delta_HR_hard_limit = -1
            b% s2% delta_HR_hard_limit = -1
            write(*,*) "delta_HR_hard_limit = ", b% s2% delta_HR_hard_limit
         end if

Cheers, Daniel

[Debraheem]

Hey Daniel, thanks, I see now where you relaxed the hard limits in the binary_extras.

I adjusted the net, but the model still does not converge. I don't want to modify your inlists away from your working example. Would you be open to resharing a working clean model directory?

on "delta_HR_limit":

MESA defaults with delta_HR_limit = -1 delta_HR_hard_limit = -1 You have: delta_HR_limit = 0.0025d0 delta_HR_hard_limit = 0.02d0 I would consider relaxing these limits, or at least relaxing the "hard_limit" as this will cause a timestep dive. Perhaps, try: delta_HR_limit = 0.0025d0 delta_HR_hard_limit = -1 This will try to make the solver limit the movement in the HR without preventing it from taking steps that exceed this limit. Movement in the HR is not always as "pretty" as published results would have everyone believe, so I wouldn't hammer too hard on HR timestep controls as there is some natural variation that can happen in the HR, especially when resolving dynamical features. Hope this helps :).

[orlox]

Hey there,

@Debraheem , I had discussed this in person with @Daniel-Pauli , and if I understood him correctly then there is a detail missing in the description of the issue. Please correct me if I'm wrong @Daniel-Pauli. I understood that you were seeing that the simulation without the hard limit runs without retries while producing no delta HR jump above the hard limit used in the other simulation. This would imply that adding the hard limit should make no difference in the two simulations.

If that is the case, then it would indicate the hard limit is activated spuriously, and in a MWE you could specify the exact time step where that happens.

[Daniel-Pauli]

Heyho,

@orlox, yeah this is what I thought in the beginning, but I did the calculations and the limit is triggered correctly.

However, I still posted it because I did not understood why the retry is working properly most of the time, but when the secondary is rotating close to critical rotation rate it leads to too small time steps and weird jumps in the HRD but when they are turned off, it just continues its way in the HRD without major issues. In an ideal world even with retrys I should get comparable solutions and not one crashing and one not. But as I understand it now, the solver is more strict and hence the behavior is different and my inlists might need an update.

I have tried now the "convergence_ignore_equL_residuals = .true." and it helps a bit, the model does not crash, but the timestep is still very low (1e-2-1e-3) years and hence it would take like forever.

As @Debraheem said, without hard limits I still can get through the mass transfer phase and I do not see major differences to my old models which were limiting the movement in the HRD more, so that is good news. However I thought I still make you aware of this issue, that the solver seems to struggle when I have an rapidly rotating accretor.

@Debraheem: If you are still interested, I can install the latest version of MESA on another machine and without any modifications and to prepare a new clean working folder so that your models will converge and you can see in more detail what is happening. But I am busy today, so I will probably do it tomorrow afternoon or so.

AHH btw. "Movement in the HR is not always as "pretty" as published results would have everyone believe". In my old version 10398 (which had a "working" hard limit, the track was prettier) thats why I was aiming for it, but I guess with the current settings my models yield results which are good enought for my purposes. Thanks for all your help and patience :)

[orlox]

One thing I can suggest that has changed since 10398 is how mass loss is treated in the energy equation (see section 3.3 of Paxton+ 2019). I have a bit of an informal feeling that one of the terms in that method can lead to jumps in surface properties, which can be bad for anything undergoing mass transfer. But I have not tested this out in detail, so take it with a grain of salt.

What I mean is the leakage term described by Eq. (26) of Paxton+ 2019, which is controlled by eps_mdot_leak_frac_factor in controls. You could try removing this with the option

eps_mdot_leak_frac_factor = 0d0

in the controls section for both stars. Doing this will affect the physics of your simulation, but if I remember correctly should provide something that is physically consistent with what was done previously on 10398. Perhaps you can try again your simulation with this change?

[Daniel-Pauli]

Hi everyone,

I have installed the latest version of MESA (r23.05.1) with the latest MESA_SDK (23.7.3) on a new machine. I reran the examples I send in my first post (with the normal nets). The models without hard limits look and behave exactly the same. However, with this installation I am also no longer able to reproduce the crashing models having the hard limits (the timestep gets small but it continues). I noticed that the MESA_SDK was released just recently and that my installation on the laptop (the one producing the crashes) has an older version (22.6.1). Now I would like to check it that is really the reason why the code crashes.. But one question: If I install the newer version of the SDK do I need to "reinstall" MESA or is it sufficient to just make another link to the new SDK in my environment?

Cheers, Daniel

[warrickball]

You should reinstall MESA with ./install but should be able to leave the unpacked data in place by running ./touch instead of ./clean. That is, presuming you've set your environment variables correctly and loaded the correct SDK, you should re-install with

$ cd $MESA_DIR
$ ./touch
$ ./install

[Debraheem]

Hello @Daniel-Pauli, were you able to resolve your issue?

[Daniel-Pauli]

Hi everyone,

sorry for the long radiosilence, I was on holidays and did not check my model calculations. Indeed the issue seems to be resolved with the new MESA SDK version. Still interesting, as the version I was using before was the latest when the MESA version r23.05.1 was released.

Thanks for all your help, Daniel

[Debraheem]

No worries, and nice to hear your issue was resolved! That is interesting indeed. Feel free to open another thread if your run into another issue.