`fe_core_infall_limit` not behaving as expected?

Describe the bug

Trying to run models to the onset of core-collapse, I encounter (unsurprisingly) lots of numerical difficulties. One particularly annoying one is that models reach the termination condition fe_core_infall_velocity not because of a large infall velocity within the Fe core, but instead because of velocity "spikes" in the outer layers (for a stripped model with no H envelope, but I've seen this in H-rich pre-SN models as well).

For example, using the trace_history_value_names to get terminal output close to core-collapse, I can find models that successfully reach the termination condition according to MESA and report:

                                 Fe_core        1.8818177594246457
                          fe_core_infall    1.1951597418765050D+04
                      non_fe_core_infall        0.0000000000000000
                               rel_E_err    1.2628869479237786D-05

and correspondingly:

stop because fe_core_infall > fe_core_infall_limit
    0.1195159742E+10    0.1000000000E+09

However, the attached pgplot shows the fe56 mass fraction profile (in cyan, right y-axis, ending at 2Msun, or, according to the terminal output 1.88 -- seems good), and the velocity profile (in yellow, left y-axis units cm/s). Clearly ,all the fe-rich material is still with v=0, but the surface (where numerical artifacts create a non-zero, and at this timestep negative, v) hits the threshold for collapse and the model terminates:

The relevant inlist lines in my setup are:

non_fe_core_infall_limit                 = 1d99
non_fe_core_infall_mass               = 1d99

fe_core_infall_limit                         = 1d8
fe_core_infall_mass                       = 1.0d0 !same behavior with default 0.1d0 or 2.0d0 (while m_fe~1.88 ?!)
when_to_stop_rtol                          = 1d-2
when_to_stop_atol                         = 1d-3

where I set non_fe_core_infall_mass and non_fe_core_infall_limit so that the model would never stop because of collapse outside of the fe core, but this does not seem to achieve the intended purpose.

Grepping the code for fe_core_infall = I see that this variable seems to be set by $MESA_DIR/star/private/report.f90, and printing (from run_star_extras.f90) the value of s% non_fe_core_infall_mass I do get 1d99 as expected. This should prevent the assignment of a value to s% non_ge_core_infall at line 426 in report.f90 because the if statement just above should be checking that mass_sum can never exceed 1d99 and indeed the terminal reports a non_fe_core_infall = 0, so far so good . However, at this point I get confused -- I suspect the problem is in this file but can't quite confirm it.

In particular: line 407 and 422 are adding to the counter mass_sum the mass coordinate s%m(k) instead of the mass within cell k (s%dm(k), note the d) which would seem more appropriate to tally the amount of mass that is collapsing by adding each cell with sufficiently high infall velocity.

Moreover, line 411 sets s% fe_core_infall to minus the minimum of the velocity (at index k_min), but it does not check whether that minimum is within the fe_core_mass mass coordinate, which seems to be the main problem here.

I'm presently experimenting with setting fe_core_infall to 1d99 and re-implement this in extras_finish_step which should be possible, but wanted to bring this up since it may affect others.

To Reproduce I have a rather complex setup (trying to make fast rotating chemically homogeneous stars), which may be why I am encountering this issue that was not noticed before. Specifically my model is by now He/C/O-rich at the surface, and he_core_mass is the total mass of the star (note that line 415 in report.f90 sets non_fe_core_mass to he_core_mass). I can send a model at Si depletion or a photo even later and my inlists/run_star_extras.f90 if useful, but none of the many complications I've introduced relate to this issue here. I did check that the core collapse test cases behave as intended.

Expected behavior If MESA states stop because fe_core_infall > fe_core_infall_limit I expect that the minimum of the velocity is lower than fe_core_infall_limit within the Fe-rich material. If I set non_fe_core_infall_limit = 1d99 I expect the model to continue no matter what the velocity profile does outside of the Fe core (or possibly crash).

Note that of course the elephant in the room is "why do the outer-layers go crazy?" I am aware/and have opinions on using some "friction" to prevent this which are being re-implemented (by @Debraheem?), but should we feed back the kinetic energy of these velocity spikes into the structure? I've also tried to implement something like velocity_q_upper_limit but in tau instead, but that is not as easy to hack from run_star_extras.f90 (I've done it finding the q value corresponding to tau that I want and updating velocity_q_upper_limit, but even setting a fairly large tau > 1d6 I couldn't manage to confine the problems so they don't impact the photosphere). Another potential problem is that the stopping condition is reached during a "retry", when the fe_core_infall goes from >0 to < -1d8 in one timestep, which should not be possible...

Desktop (please complete the following information):

I doubt this is relevant

OS: Debian 6.1.67-1 (2023-12-12) x86_64 GNU/Linux
Version: r23.05.1

For completeness, here is a snippet of code in my run_star_extras.f90 (inextras_finish_step) that attempts to reimplementreport.f90` -- unfortunately my model then crashes (because of the surface going crazy) before the Fe core starts significantly infalling, so as of now this is still untested:

    if (s% fe_core_mass > 0) then
       m_infall = 0.0d0
       ! find index outer edge Fe core
       s% xtra6_array = s%m/Msun - s%fe_core_mass
       k_fe = minloc(abs(s% xtra6_array(1:s%nz)), dim=1)
       do k=k_fe, s%nz !loop from outer edge inward
          if (-s%v(k) >= 1d8) m_infall = m_infall + s%dm(k)
       end do
       if (m_infall >= 0.1d0 * Msun) extras_finish_step = terminate
    end if

(which assumes v_flag=.true. and only considers the velocity and mass inside the fe_core_mass set elsewhere)

I confirm that I can obtain the velocity profile I expect by setting fe_core_infall_velocity = 1d99 in the inlist (i.e., turn that off completely) and re-implementing it in a fashion similar to above:

(note: to not have this particular model fail before reaching the desired termination condition the loop above is not only inside the Fe core, but inside the outermost layer for which the sound crossing time to the edge of the Fe core is smaller or equal to the current timestep)

I could open a PR with an updated report.f90, where I think the minimum change necessary is line 410-412 which should check if k_min >= k(edge of Fe core), but I'll wait for the opinion of the developers.

@mathren , thanks for documenting this! It does seem like a bug in report.90, and it seems you've confirmed your suspicions using a run_star_extras. If you'd like to open a pr, that'd be wonderful! I will merge it into a branch to test.

Note that of course the elephant in the room is "why do the outer-layers go crazy?" I am aware/and have opinions on using some "friction" to prevent this which are being re-implemented (by @Debraheem?), but should we feed back the kinetic energy of these velocity spikes into the structure?

The most recent release candidate r24.02.1-rc1 includes the velocity damping from r15140, although the drag_energy was off by default. In the candidate release, however, we do incorporate the drag energy back into the energy equation

295: sources_ad = eps_nuc_ad - non_nuc_neu_ad + extra_heat_ad + Eq_ad + RTI_diffusion_ad + drag_energy

I think it would be useful to have a control that toggles drag_energy on or off. We could incorporate this into your pr.

Happy yo give it a try! I may need some pointers to implement the toggle for the drag_energy (I know I've played with this but failed in the past). While we are at it, what do you think about velocity_tau_upper_bound? Picking a certain q that would work for both RSG and stripped progenitors can be tricky (in my experience), but velocity_q_upper_bound is acting quite deep and it's not easy to hack form run_star_extras.f90...

Happy yo give it a try! I may need some pointers to implement the toggle for the drag_energy (I know I've played with this but failed in the past)

No worries, I can help at any point in the process.

While we are at it, what do you think about velocity_tau_upper_bound? Picking a certain q that would work for both RSG and stripped progenitors can be tricky (in my experience), but velocity_q_upper_bound is acting quite deep and it's not easy to hack form run_star_extras.f90...

If I understand you correctly, you'd like to implement a control similar to velocity_q_upper_bound as velocity_tau_upper_bound? While I think this could be done, I don't know if it will solve your issue. In my experience turning the velocity off entirely in the envelope has resulted in even stranger things. Nonetheless, certainly feel free to try it out!

MESAHub / mesa

`fe_core_infall_limit` not behaving as expected? #626