remove/modify fe_core_mass limit for tracking velocity

[Debraheem]

Currently, during Fe core infall the peak velocity will form off center during the initial stages of CC and then gradually move inward into the Fe core as the infall progresses. If we are trying to track the infall speed as CC proceeds, we must be lenient in our definition of core infall. The current infall condition only considers the Fe core to be infalling if the peak velocity of the infall is inside the Fe core. This results in discontinuous jumps in the Fe_core_infall condition, making it sensitive to the infall profile.

One simple remedy is to remove the Fe-core infall limit and revert this line of code
(s%m(k_min) <= s%fe_core_mass*msun)) then
from #628.

Another solution would be to relax the condition to use the si_core_mass instead, however we don't currently track the si_core_mass by default in MESA any more.

[mathren]

The reason I added that line was to prevent models with infall velocity larger than fe_core_infall_limit only outside the fe core from triggering the termination condition, which was happening (sometimes not even due to infall of the core, but just numerically spurious velocities). In these situation the fe_core_infall would be set to velocity(k_min) even if k_min is outside the fe_core which is incorrect. I think we should still prevent those situations, because they are easy to miss unless one inspects the velocity profile closely.

Typically the "infall" condition should be just a few hundred millisec from core-bounce for the models to be useful as CCSN progenitors. MESA of course cannot go all the way to core bounce because of missing physics (in the neutrino interaction and EOS), so it can't "predict" when core-bounce will be and stop a bit earlier. The condition fe_core_infall_limit <= 1000km/s from Woosley et al. 2002 is designed to be close enough to core-bounce.

In my opinion, if the velocity profile in the Fe core is not infalling fast enough, what happens outside the Fe core should not stop the model. On the other hand, maybe demanding that the minimum of the velocity profile is in the Fe core is too restrictive, a condition could be "the minimum velocity within the fe core should be below the threshold" (regardless of whether it is the absolute minimum of the velocity throughout the star or not).

In a sense, a few lines above your edit: k_min = minloc(velocity(1:nz), dim=1) should probably only check the minimum within the fe_core instead of the whole array (1:nz) (which would also make the if statement you removed useless at that point). This way fe_core_infall conditions only depend on variables within the fe_core (and are not checked until there is a fe_core), and nothing outside it?

[Debraheem]

In a sense, a few lines above your edit: k_min = minloc(velocity(1:nz), dim=1) should probably only check the minimum within the fe_core instead of the whole array (1:nz) (which would also make the if statement you removed useless at that point). This way fe_core_infall conditions only depend on variables within the fe_core (and are not checked until there is a fe_core), and nothing outside it?

I agree. We can check for the max velocity in the fe-core instead of asking for the max_velocity and then checking that the max velocity is inside the fe-core.

[mathren]

To this end, I've used this snippet in extras_finish_step in my run_star_extras.f90:

    ! ! custom Fe core collapse
    if (s% fe_core_mass >0.0d0) then
       !log more stuff in the terminal
       s% num_trace_history_values     = 5
       s% trace_history_value_name(1) = 'Fe_core'
       s% trace_history_value_name(2) = 'fe_core_infall'
       s% trace_history_value_name(3) = 'non_fe_core_infall'
       s% trace_history_value_name(4) = 'rel_E_err'
       s% trace_history_value_name(5) = 'log_rel_run_E_err'
       s% trace_history_value_name(6) = 'dt_div_max_tau_conv'
       k = s%nz
       do while (s% m(k) <= s% fe_core_mass * Msun)
          k = k-1 ! loop outwards
       end do
       ! k is now the outer index of the fe core
       if (maxval(abs(s%v(k:s%nz))) >= s% fe_core_infall_limit) then
          s% termination_code = t_fe_core_infall_limit
          write(*, '(/,a,/, 99e20.10)') &
               'stop because fe_core_infall > fe_core_infall_limit', &
               s% fe_core_infall, s% fe_core_infall_limit
          print *, "treshold v used", maxval(abs(s%v(k:s%nz)))
          extras_finish_step = terminate
       end if
    end if

[mathren]

@Debraheem should we use instead of

s% fe_core_infall = - minval(s%v(k:nz))

the following:

s% fe_core_infall = abs(minval(s%v(k:nz))

so that fe_core_infall is always defined positive (which we want), and if the minimum is positive for some reason we don't get a negative number?

[Debraheem]

If the minimum is positive, the fe_core_infall stopping condition won't be toggled, because s% fe_core_infall = - minval(s%v(k:nz)) will ensure s% fe_core_infall is negative.

the infall condiiton is toggled when s% fe_core_infall > s% fe_core_infall_limit, where fe_core_infall_limit = (+)

if you choose s% fe_core_infall = abs(minval(s%v(k:nz)), you risk prematurely exiting from a positive velocity, like from an outward moving shock.

I prefer the former.