Feature/superad_reduction: opt-in turnover-time limiter with v_c floor and smooth relaxation

[matteocantiello]

Summary

This PR adds an opt-in physical refinement to the existing
use_superad_reduction machinery (Jermyn+23), making the radiative-gradient
throttle aware of both the local convective response time and its
previous-step state.

The PR introduces one new feature, together with two coupled improvements that
were identified during validation:

1. a turnover-time limiter,
2. a velocity floor for slow-convection cells, and
3. a smooth-relaxation formula anchored to the previous converged throttle.

All new behavior is opt-in. With defaults, the existing test_suite output is
bit-for-bit unchanged.

What this PR does

1. Turnover-time limiter

New control:

superad_reduction_use_turnover_limit

Type/default:

logical, default .false.

When enabled, the throttle excess

Gamma_fac - 1

is multiplied by

f_tau = 1 - exp(-dt / tau_conv)

where

tau_conv = scale_height(k) / mlt_vc_old(k)

This leaves the throttle essentially unchanged when

dt >> tau_conv

but smoothly attenuates it when

dt < tau_conv

that is, when convection cannot respond as quickly as the implicit step is
asking it to.

2. Velocity floor for slow-convection cells

New companion control:

superad_reduction_turnover_vc_floor_frac

Type/default:

real, default 0d0

Before computing tau_conv, this floors mlt_vc_old(k) at

frac * csound_face(k)

This is needed in radiation-pressure-dominated layers, especially in the
iron-bump surface zone, where MLT can significantly under-predict the relevant
response speed. Even when MLT predicts very small convective velocities,
perturbations can still communicate through the layer on an acoustic timescale.

Without this floor, mlt_vc_old(k) can become so small that

f_tau -> 0

which fully suppresses the throttle in exactly the cells where radiation-pressure
inversions make the Newton solve fragile.

3. Smooth-relaxation formula

The original limiter formula was

Gamma_fac_new = 1 + f_tau * (Gamma_fac_inst - 1)

This relaxes Gamma_fac toward 1, i.e. toward no throttle, when
dt << tau_conv.

That is not the desired behavior when the throttle is already active. In that
case, the throttle should also not be released faster than convection can
respond. Abruptly relaxing toward Gamma_fac = 1 can produce a discontinuous
structural rearrangement, visible as a small HRD jump near the end of the
to_cc phase when the timestep collapses during iron-core formation.

This PR replaces that expression with a first-order relaxation anchored to the
previous step's converged value:

Gamma_fac_new = Gamma_fac_old + f_tau * (Gamma_fac_inst - Gamma_fac_old)

or equivalently,

Gamma_fac_new = (1 - f_tau) * Gamma_fac_old + f_tau * Gamma_fac_inst

This is the discrete form of

dGamma_fac/dt = (Gamma_fac_inst - Gamma_fac) / tau_conv

The limiting behavior is then:

• dt >> tau_conv: f_tau -> 1, so Gamma_fac -> Gamma_fac_inst
• dt << tau_conv: f_tau -> 0, so Gamma_fac -> Gamma_fac_old

No additional control is needed. On the first step, when Gamma_fac_old is not
yet available, the guard falls through and the limiter reduces to the original
first-step behavior. On subsequent steps, the smooth-relaxation form is used.

Implementation

Main implementation changes:

File	Change
star_data/public/star_data_step_work.inc	Adds superad_reduction_factor_old(:) and have_superad_reduction_factor, following the existing paired _old array pattern.
star_data/private/star_controls.inc	Declares the two new namelist controls.
star/private/ctrls_io.f90	Wires the controls into &controls, including namelist, store, and write support.
star/defaults/controls.defaults	Adds defaults and documentation. Both controls default to off.
star/private/alloc.f90	Allocates superad_reduction_factor_old in star_info_old_arrays, avoiding null-pointer issues in check_sizes.
star/private/evolve_support.f90	Copies superad_reduction_factor into superad_reduction_factor_old at the end of each successful step, mirroring the mlt_vc -> mlt_vc_old pattern.
star/private/init.f90 and star/private/remove_shells.f90	Initialize/reset the have_superad_reduction_factor flag.
star/private/turb_support.f90	Adds the limiter block inside set_superad_reduction, including the velocity floor, smooth-relaxation formula, and guard logic.

Validation

A. 60 Msun A/B/C benchmark

Restart from a photo shortly before the He-exhaustion timestep collapse, then
run through xa_central_lower_limit.

Run	End condition	Steps	Retries
limiter off, baseline	xa_central_lower_limit	1940	14
limiter on, no floor	max_model_number	6000	780
limiter on, floor = 1d-3 * cs	xa_central_lower_limit	1910	15

The bare limiter, without the velocity floor, suppresses the throttle in roughly
5-17 surface iron-bump cells. Radiation-pressure inversions then reappear,
leading to many retries and failure to reach central He exhaustion.

With the 1d-3 * cs floor, those cells remain regularized and the run returns
to baseline-quality numerics while keeping the physical limiter active.

B. Smooth relaxation removes the late-phase HRD jump

For the 60 Msun to_cc phase, restarted from after_core_c_burn:

Formula	Endpoint log T_eff	Endpoint log L	fe_core_infall_limit?
old: Gamma_fac = 1 + f_tau * (Gamma_fac_inst - 1)	3.58	5.58	yes
new: Gamma_fac = Gamma_fac_old + f_tau * (Gamma_fac_inst - Gamma_fac_old)	3.70	6.11	yes

Both runs reach fe_core_infall_limit at the same central T_c / rho_c
state. The smooth-relaxation formula keeps the surface evolution continuous to
the endpoint and removes the visible HRD jump.

C. Full pipelines to core collapse at solar metallicity

With the full feature stack enabled:

superad_reduction_use_turnover_limit = .true.
superad_reduction_turnover_vc_floor_frac = 1d-3

and using the smooth-relaxation formula:

Mass	Outcome
30 Msun	fe_core_infall_limit reached
40 Msun	fe_core_infall_limit reached
60 Msun	fe_core_infall_limit reached

Matched controls without the limiter stalled during He burning or earlier.

Backwards compatibility

Both new controls default to off.

With

superad_reduction_use_turnover_limit = .false.

the new limiter code path is skipped, and test_suite output is bit-for-bit
identical to the parent branch.

On the first step with the limiter enabled, have_superad_reduction_factor is
false, so the smooth-relaxation formula falls through and the limiter produces
the same first-step output as the original two-term form. From the second step
onward, the smooth-relaxation form is active.

[Debraheem]

Amazing! I think we should refactor this a bit but this is great. I am however a little unconvinced that we should be using an exponential to limit the enhancement. This is one of those cases where I think linear might be more self consistent, and stable, but maybe I can be proven wrong.

[aurimontem]

Adding
superad_reduction_use_turnover_limit = .true.
superad_reduction_turnover_vc_floor_frac = 1d-3

works great (no crash, no excursions, no spurious oscillations which arise with MLT++ and superad in hydro mode without the limiter) for our hackathon test case with O-burning 25Msun that crashes with excursions w/ and w/o hydro.

Excellent stuff!!!

[matteocantiello]

Latest updates

This PR now includes three additional fixes needed to make the smooth
superad-reduction limiter robust across mesh changes, retries, and throttle
release.

1. Mesh remapping of the previous-step superad factor

The previous-step superad-reduction factor is now carried through mesh remaps,
mirroring the existing mlt_vc / mlt_vc_old handling.

This adds the prev_mesh_superad_reduction_factor array and the associated
plumbing through:

• alloc.f90
• evolve.f90
• mesh_adjust.f90
• adjust_mesh.f90
• remove_shells.f90

The remap path includes the needed save/restore, interpolation, pruning, and
caller updates. This is required because the smooth limiter uses the previous
step's converged superad-reduction factor; after a mesh change, that quantity
must remain defined on the current mesh.

2. Correct throttle-release behavior

The limiter block has been moved outside the

if (Gamma_term > 0d0)

block.

The guard has also been changed from requiring

Gamma_factor > 1d0

to allowing the limiter to operate when either the new step or the previous step
is throttled:

Gamma_factor > 1d0 .or. Gamma_factor_old > 1d0

This matters because the limiter must control both directions:

• turning the throttle on when the current step needs it, and
• releasing the throttle only on the relevant response timescale when the
  previous step was throttled.

Without this, the code could still release the throttle too abruptly when the
instantaneous new value no longer required superad reduction.

3. Snapshot protection

This PR now adds an okay_to_set_superad_reduction_factor flag, mirroring the
existing okay_to_set_mlt_vc pattern.

The flag is set to false during:

• initialization,
• start of step,
• retries,
• shell removal,

and is set to true only after new_generation snapshots are valid.

The write at the end of set_superad_reduction is gated on this flag.

This prevents invalid or intermediate values of the superad-reduction factor
from being saved and later used as the previous-step reference value. In other
words, the smooth limiter now only relaxes against a clean, converged snapshot.

[Debraheem]

we might want to mimic the behavior of public star_utils tau_conv functions for this:

mesa/star/private/star_utils.f90

Lines 3698 to 3751 in aa27a08

	real(dp) function conv_time_scale(s,k_in) result(tau_conv)
	type (star_info), pointer :: s
	integer, intent(in) :: k_in
	integer :: k
	real(dp) :: brunt_B, alfa, beta, rho_face, Peos_face, chiT_face, chiRho_face, &
	f, dlnP, dlnT, grada_face, gradT_actual, brunt_N2
	if (.not. s% calculate_Brunt_B) then
	tau_conv = 0d0
	return
	end if
	k = max(2,k_in)
	brunt_B = s% brunt_B(k)
	call get_face_weights(s, k, alfa, beta)
	rho_face = alfa*s% rho(k) + beta*s% rho(k-1)
	Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1)
	chiT_face = alfa*s% chiT(k) + beta*s% chiT(k-1)
	chiRho_face = alfa*s% chiRho(k) + beta*s% chiRho(k-1)
	f = pow2(s% grav(k))*rho_face/Peos_face*chiT_face/chiRho_face
	dlnP = s% lnPeos(k-1) - s% lnPeos(k)
	dlnT = s% lnT(k-1) - s% lnT(k)
	grada_face = alfa*s% grada(k) + beta*s% grada(k-1)
	gradT_actual = safe_div_val(s, dlnT, dlnP) ! mlt has not been called yet when doing this
	brunt_N2 = f*(brunt_B - (gradT_actual - grada_face))
	if(abs(brunt_B) > 0d0) then
	tau_conv = 1d0/sqrt(abs(brunt_N2))
	else
	tau_conv = 0d0
	end if
	end function conv_time_scale
	subroutine set_conv_time_scales(s)
	type (star_info), pointer :: s
	integer :: k
	real(dp) :: tau_conv
	include 'formats'
	s% min_conv_time_scale = 1d99
	s% max_conv_time_scale = 0d0
	do k=1,s%nz
	if (s% X(k) > s% max_X_for_conv_timescale) cycle
	if (s% X(k) < s% min_X_for_conv_timescale) cycle
	if (s% q(k) > s% max_q_for_conv_timescale) cycle
	if (s% q(k) < s% min_q_for_conv_timescale) exit
	tau_conv = conv_time_scale(s,k)
	if (tau_conv < s% min_conv_time_scale) &
	s% min_conv_time_scale = tau_conv
	if (tau_conv > s% max_conv_time_scale) &
	s% max_conv_time_scale = tau_conv
	end do
	if (s% max_conv_time_scale == 0d0) s% max_conv_time_scale = 1d99
	if (s% min_conv_time_scale == 1d99) s% min_conv_time_scale = 0d0
	end subroutine set_conv_time_scales

[Debraheem]

My previous comment and the following comment only apply if we decide we'd like to make this an explicit quantity that is fixed over solver iterations. Right now you're mixing an implicit scale height with an explicit start of step velocity, and so tau_conv is being recalculated every solver iteration and changing. We should think about if that is desirable or if tau_conv should just be a fixed value over solver iterations.

if we go the explicit path:
We can just compute this quantity once at the start of step so it doesn't need recalculated. Then the simple way to do this is probably to just add a star pointer variable s% tau_conv and set s% tau_conv(k) = tau_conv inside of set_conv_time_scales which is already called with brunt at the start of step here

mesa/star/private/hydro_vars.f90

Lines 545 to 553 in aa27a08

	if (.not. skip_grads) then
	if (dbg) write(*,*) 'call do_brunt_B'
	call do_brunt_B(s, nzlo, nzhi, ierr) ! for unsmoothed_brunt_B
	if (failed('do_brunt_B')) return
	if (dbg) write(*,*) 'call set_grads'
	call set_grads(s, ierr)
	if (failed('set_grads')) return
	call set_conv_time_scales(s) ! uses brunt_B
	end if

Then you can just check for s%dt/s%tau_conv(k) here inside turb_support.

[Debraheem]

If convergence is unaffected, the explicit path should be faster, and avoids recalculating things that don't change over the iterations.

[Debraheem]

If the zone begins the step as radiative, the limiter is not applying. This needs addressed somehow. I think the best solution is to use a different timescale that doesn't depend on mlt_vc at all. Wheteher it's the 1/N^2 or a sound_crossing/dynamical time.

[Debraheem]

if we think of a better timescale, we won't need an axuillary control for flooring V. Either way, we might want the floor to a parameter rather than a control, since it's not really something that should be changing, rather it's just clipping bad values.

[Debraheem]

these comments are superfluous here. Let's keep the docs in the contols, and keep a brief, concise but clear explanation of what is happening in the superadreduction scheme when it's called.

[Debraheem]

also superfluous comments that should be removed.

[Debraheem]

Is this routine ever called with dt = 0? does mesa ever call mlt dt = 0? maybe during the pre-ms model builder?

[Debraheem]

I don't understand this statement "! Clamp at >= 1 so we never relax from a non-physical sub-1 anchor."

Shouldn't max(s% superad_reduction_factor_old(k), 1d0) actually be min(s% superad_reduction_factor_old(k), 1d0)? the factor should always be between 0 and 1. And remeshing wouldn't make superad_reduction_factor_old > 1, if it's never > 1 on the previous mesh.

[Debraheem]

perhaps my confusion here is that we are using gamma_factor instead of grad_scale. I think we should be limiting gradscale not gamma_factor. gradscale is always between 0 and 1

grad_scale = (gradr-gradL)/(Gamma_factor*gradr) + gradL/gradr

[Debraheem]

gamma factor can be large and positive

[Debraheem]

Or more specifically i think we should be limiting 1/gamma_factor.

[Debraheem]

okay i think this is just a minor bug, here i explain below.

The code uses:

gradr_scaled = gradL + (gradr - gradL)/Gamma_factor

or equivalently:

gradr_scaled = grad_scale * gradr

where:

grad_scale = gradL/gradr + (1/Gamma_factor)*(1 - gradL/gradr)

So Gamma_factor is not the quantity that multiplies gradr going into MLT/TDC, rather Gamma enteres grad_scale in the denominator:

eta = 1/Gamma_factor

because these multiply the effective radiative gradient input.

Since our goal is to limit changes in gradr on a tunrover time, I think we should be limiting eta = 1/gamma_factor or grad_scale, not Gamma_factor itself.

Example:

Gamma_old = 1
Gamma_inst = 100 ! instantaneous gamma that superad wants with no limiter.
f_tau = 0.1

Current PR limiter:

Gamma_new = Gamma_old + f_tau*(Gamma_inst - Gamma_old)
          = 1 + 0.1*(100 - 1)
          = 10.9

eta_new = 1/Gamma_new = 0.092

That means about 91% of the reduction is applied immediately, even though f_tau is only 0.1.

If instead we limit the effective gradr scaling:

eta_old = 1/Gamma_old = 1
eta_inst = 1/Gamma_inst = 0.01

eta_new = eta_old + f_tau*(eta_inst - eta_old)
        = 1 + 0.1*(0.01 - 1)
        = 0.901

Gamma_new = 1/eta_new = 1.11

That applies only about 10% of the reduction, which matches the meaning of f_tau.

The easy soution is just to use use the limiter on eta = 1/Gamma_factor, e.g. using eta:

eta_inst = 1/Gamma_factor ! current iterate gamma_factor before limiting. 
eta_old  = 1/Gamma_factor_old

[aurimontem]

So by way of an update on the adapted MESA VI S7.2 25M to O depletion model, without the limiter the model stalls whether or not hydro is on, and without hydro it’s worse. With MLT++ it does not stall but it does produce ratty behavior. With no superad/mlt++, the evolution itself is OK, and it gets to O depletion, though the physicality of MLT itself on the RSG branch is … dubious. With the initial version of the limiter, it jumps (almost instantaneously) from the superad HR track to the standard MLT track and the rattiness picks up at the end — BUT it gets to O depletion fairly easily. With the new version of the limiter where you don’t turn the reduction off, but rather keep it fixed for shorter dt, it gets to O depletion fairly easily, but the evolution in the HR gets very ratty.

Might push to core-collapse to diagnose further (and see how much the core evolution is affected)

[Debraheem]

So by way of an update on the adapted MESA VI S7.2 25M to O depletion model, without the limiter the model stalls whether or not hydro is on, and without hydro it’s worse. With MLT++ it does not stall but it does produce ratty behavior. With no superad/mlt++, the evolution itself is OK, and it gets to O depletion, though the physicality of MLT itself on the RSG branch is … dubious. With the initial version of the limiter, it jumps (almost instantaneously) from the superad HR track to the standard MLT track and the rattiness picks up at the end — BUT it gets to O depletion fairly easily. With the new version of the limiter where you don’t turn the reduction off, but rather keep it fixed for shorter dt, it gets to O depletion fairly easily, but the evolution in the HR gets very ratty.

Might push to core-collapse to diagnose further (and see how much the _core_ evolution is affected)

maybe some of the fixes i suggested will help remove this ratty hr behavior.