Problem with eps_mdot and high mass loss rates

[jschwab]

@orlox, when you escape grant writing, please provide an example of the case that drove you to resurrect some of the old eps_grav code in your run star extras.

[jschwab]

@orlox Do you have any example of this issue?

[orlox]

files.zip Here's one example, includes saved model and inlists. Specific model here is a 30 Msun star at low metallicity that I strip at a high rate of 1 Msun/yr. Files work with 0292d48, after about 200 steps where the model advances steadily, outermost layers go bonkers. In 15140 I can circumvent this by using the old eps_grav homologous/non-homologous form, with which the model goes forward without any issues. Can also provide an example of this in 15140 with the old eps_grav, more recent versions don't have dlnd_dt_const_q and dlnd_dt so it's a bit more cumbersome to implement.

PS: grant writing has eaten way more of my time than I expected, but I'll finally be free after this week.

[adamjermyn]

Thanks!

This behavior is caused by the 'leak' corrections in eps_mdot. Briefly, when the thermal time in a region becomes longer than the mass-loss timescale (xm/mdot) eps_mdot accounts for the fact that material does not adjust to the local temperature profile as it moves upwards through the star and instead retains its starting entropy.

You can see this by setting

    eps_mdot_leak_frac_factor = 0d0

in the controls section of your inlist. This turns off the correction and tells eps_mdot 'I know this is physically wrong but I really just want to get a model running'. The model then happily runs through at least model 3200 without any of that noise in the surface temperature profile. It does eventually get pushed down to very short dt ~ 1e-6yr around model 3300 (at which point the mass is down to 12.5 Msun, the core is mostly Oxygen, and it's doing an awful lot of burning). Not sure if that's acceptable, but I'm guessing that's similar to what you were seeing with eps_grav?

If you add an entropy panel to the plots you can see that near/in/above the non-isentropic regions are where this leakage correction makes a difference in eps_mdot, and that's exactly what you expect if you're moving material faster than the local thermal time.

More generally, the leakage correction often makes models harder to run. I suspect there are a few reasons for this:

1. The actual physics is harder to model! eps_grav got around this by being pretty wrong in this limit. eps_mdot is much closer to right in this limit, but that sometimes means the problem is harder and often it means you need to take smaller timesteps.
2. When Mdot < 0, mass is moving from a region where MESA generally uses lower mass resolution to a region where it needs higher resolution. It's not too surprising that a low-resolution entropy profile superimposed on an atmosphere that 'wants' higher resolution runs into trouble.
3. The way we achieve this adiabatic motion is sort of funny. We do it by letting the material change state to the local (rho, T) profile, then we give it the amount of heat it needs to get back to its starting entropy. If the solver did a perfect job and the equations were perfect and numerics were perfect this would probably work better than it does. I'm not sure how much this contributes, but I wouldn't be surprised if this mattered...

Anyway, I don't see anything to suggest a problem in eps_mdot here. If you're in a position where you just need to run a model and don't actually care about getting the energetics right I'd say just turn the leak factor to zero and run it. If you do actually care about the energetics you're in a bit of trouble here because it looks like a hard problem that eps_grav wasn't appropriate for, and with eps_mdot it'll take a lot of handholding to get it to converge. In particular I think this limit will force you to use a timestep comparable to (xm/mdot) in the shallowest radiative zone, which is often going to be very small. But I don't know of a better way to implement this physics in a fully Lagrangian code.

Why don't you give it a whirl with the leak factor set to zero and let me know if that's good enough?

[adamjermyn]

I should add: I was able to get a model that 'converged' reasonably using additional timestep controls without disabling the leak correction. The controls were:

    delta_lgT_limit = 3d-2
    delta_lgT_limit_min_lgT = 0d0

This forces the time steps to be about 10x smaller than they were originally but produces a model that's well-behaved everywhere except right near the photosphere. If you follow on an HR diagram you'll see logL and logTeff jump around a lot, but that's purely a 'surface' effect. Deeper down everything is well-behaved and the solver doesn't complain too much.

[adamjermyn]

Some of what I said above is wrong. Setting

    eps_mdot_leak_frac_factor = 0d0

is equivalent to forcing the system to be adiabatic, which should make things worse numerically. Setting

    eps_mdot_leak_frac_factor = (enormous)

is equivalent to forcing the system to be isothermal, which should recover the older prescription.

@jschwab just sent me another case where there's something funny going on, so I'll dig in more.

[adamjermyn]

Ok. Setting the leak fraction to zero just dramatically reduces eps_mdot (because material carries most of the heat with it on its way out of the model). So that's why that improved things.

The rest of what I said holds though. When I set

   eps_mdot_leak_frac_factor = 1d99

that forces it to replicate the old non-Lagrangian eps_grav term (more or less). This gives a better-behaved result (though a less physical one!):

Let me know if you'd like me to think more about this/what you'd like me to think about...

[orlox]

I want to do some more tests on my side before closing this, the info you've given me Adam is pretty useful, thanks for that!

[adamjermyn]

Ok! Please let me know if you encounter ongoing weirdness. I want to dig into eps_mdot more anyway because of some funny behavior Josiah pointed out. It could be working as expected but it could be a bug, so the more info I get the better...

[adamjermyn]

I did a bit more digging and didn't see anything wrong, so nothing more here on my end. Do let me know if you see more weirdness!