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

docs/source/changelog.rst

@@ -30,6 +30,14 @@ New Features
 
 MESA no longer stops when reactions for which special rates are set are not in the nuclear network, only a warning is printed. This is intended to make it easier to test various network sizes without having to also change the list of special reactions.
 
+Two new profile columns, ``superad_reduction_Lrad_div_Ledd`` and ``superad_reduction_trigger``, and one new history column, ``num_cells_with_superad_reduction``, expose what the ``use_superad_reduction`` throttle is doing per cell (which threshold fired and at what proximity to Eddington). The existing ``superad_reduction_factor`` profile column and ``max_superad_reduction_factor`` history column are unchanged. The new ``superad_reduction_trigger`` is an integer flag: 1 = Eddington proximity (``Lrad/Ledd > Gamma_limit``); 2 = density-inversion criterion (Joss et al. 1973, Paxton, Cantiello et al. 2013 eq. 17); 3 = both.
+
+Opt-in convective-turnover-time limiter for ``use_superad_reduction``: when the new control ``superad_reduction_use_turnover_limit`` is set to ``.true.``, the throttle excess ``(Gamma_factor - 1)`` is multiplied by ``f_tau = 1 - exp(-dt / tau_conv)``, where ``tau_conv = scale_height(k) / mlt_vc_old(k)`` is the local convective turnover time built from the previous-step converged convective velocity. The limiter recovers the current behavior when ``dt >> tau_conv`` and smoothly suppresses the throttle when ``dt << tau_conv`` (where the implicit-step throttle would be acting faster than convection itself can respond). Default ``.false.``; bit-for-bit unchanged on every shipped test_suite case.
+
+Companion control to the turnover-time limiter: ``superad_reduction_turnover_vc_floor_frac`` (real, default ``0d0``). When set to a small positive value, floors ``mlt_vc_old(k)`` at ``frac * csound_face(k)`` before computing ``tau_conv``. This prevents ``f_tau`` from collapsing to ~0 in slow-convection iron-bump cells where the bare limiter would fully suppress the throttle and re-expose the radiation-pressure inversions to the Newton solver. With ``frac = 1d-3`` the floor restores baseline-quality numerics through the post-He burning crunch (60 M_sun benchmark: 15 retries through fe_core_infall, matching the limiter-off baseline's 14 retries; the bare limiter accumulates ~780 retries in the same window). Default ``0d0`` disables the floor and recovers the bare-limiter behavior exactly.
+
+The turnover-time limiter formula now anchors the relaxation at the previous step's converged ``Gamma_factor`` rather than at 1 (no throttle): ``Gamma_factor_new = Gamma_factor_old + f_tau * (Gamma_factor_inst - Gamma_factor_old)`` instead of the earlier ``Gamma_factor_new = 1 + f_tau * (Gamma_factor_inst - 1)``. The two limits remain physically correct (``dt >> tau_conv`` recovers the instantaneous value, ``dt << tau_conv`` holds the throttle at its previous value) but the new formula no longer asks convection to *release* the throttle on timescales shorter than ``tau_conv`` -- which the old formula did, and which caused a small but visible HRD jump at the end of the to_cc phase when ``dt`` collapsed for iron-core formation. The change required a new paired ``superad_reduction_factor_old(:)`` array (mirroring the ``mlt_vc/mlt_vc_old`` pattern). On the very first step ``Gamma_factor_old`` is unset and the formula falls through to the old behavior, preserving first-step output bit-for-bit. No new control: the new formula is strictly more correct on every step after the first and is enabled whenever ``superad_reduction_use_turnover_limit = .true.``.
+
 .. _Bug Fixes main:
 
 Bug Fixes
@@ -45,6 +53,10 @@ The parameter ``report_max_infall_inside_fe_core`` was ignored in versions r25.1
 
 ``fe_core_infall_limit`` now obeys ``when_to_stop_rtol`` and ``when_to_stop_atol`` again (broken since r11532).
 
+In ``set_superad_reduction`` (``star/private/turb_support.f90``) the density-inversion contribution to ``Gamma_term`` was being scaled by ``superad_reduction_Gamma_limit_scale`` instead of the intended ``superad_reduction_Gamma_inv_scale``, making the latter control silently inert. Restored to the as-published algorithm (Jermyn et al. 2023 eq. 64, ``alpha_2`` term). Existing test_suite cases that use ``superad_reduction`` all set ``Gamma_limit_scale = Gamma_inv_scale``, so they are bit-for-bit unaffected; the fix is observable only when the two scales are set to different values.
+
+The default value of ``superad_reduction_limit`` is changed from ``-1d0`` (uncapped) to ``100d0`` (logistic saturation at ``Gamma_factor ~= 100``). The cap was already documented as the recommended setting, but the shipped default disabled it. With the cap inactive, the raw ``1/sqrt(beta)`` boost in radiation-pressure-dominated cells can drive ``Gamma_factor`` arbitrarily large in a single Newton iteration. Existing test_suite cases that exercise ``superad_reduction`` all explicitly set this control to ``-1d0``, so they are bit-for-bit unaffected.
+
 
 .. note:: Before releasing a new version of MESA, move `Changes in main` to a new section below with the version number as the title, and add a new `Changes in main` section at the top of the file (see ``changelog_template.rst``).
 

star/defaults/controls.defaults

@@ -2655,27 +2655,168 @@
 
       ! use_superad_reduction
       ! ~~~~~~~~~~~~~~~~~~~~~
+
+      ! Implicit alternative to ``okay_to_reduce_gradT_excess`` (MLT++).
+      ! When ``.true.``, throttles the radiative gradient ``gradr`` in
+      ! convective cells that either (a) carry a radiative luminosity
+      ! approaching the Eddington limit, or (b) sit in the
+      ! density-inversion regime (Joss--Salpeter--Ostriker). MLT/TDC is then
+      ! re-solved on the throttled ``gradr``, so the final ``gradT`` is
+      ! self-consistent with MLT for the reduced radiative input.
+      ! Acts on the *input* to MLT, not the output --- unlike MLT++.
+
+      ! ::
+
+      use_superad_reduction = .false.
+
+
       ! superad_reduction_Gamma_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Radiative-Eddington trigger threshold:
+      ! ``Lrad/Ledd = (4 a T^4 / 3 P) * gradT``.
+      ! When this ratio exceeds ``Gamma_limit`` in a convective cell,
+      ! the Eddington-proximity contribution to the throttle activates.
+
+      ! ::
+
+      superad_reduction_Gamma_limit = 0.5d0
+
+
       ! superad_reduction_Gamma_limit_scale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Sensitivity of the throttle to ``Lrad/Ledd`` excess above
+      ! ``Gamma_limit``. The Eddington contribution to ``Gamma_factor``
+      ! grows quadratically (then linearly) with the fractional excess;
+      ! ``Gamma_limit_scale`` is the proportionality constant. Larger
+      ! values produce stronger throttling for a given excess.
+
+      ! ::
+
+      superad_reduction_Gamma_limit_scale = 5d0
+
+
       ! superad_reduction_Gamma_inv_scale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Sensitivity of the throttle to the density-inversion criterion
+      ! ``Lrad/Ledd > 4(1-beta)/(4-3 beta)``, where ``beta = Pgas/Ptot``.
+      ! Analogous to ``Gamma_limit_scale`` but for the inversion term.
+      ! Tune independently when inversion proximity dominates over
+      ! Eddington proximity (or vice versa).
+
+      ! ::
+
+      superad_reduction_Gamma_inv_scale = 5d0
+
+
       ! superad_reduction_diff_grads_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Superadiabaticity floor below which the throttle is silenced.
+      ! When ``gradT - gradL < diff_grads_limit``, the cell is treated
+      ! as essentially adiabatic and the throttle contribution drops
+      ! smoothly to zero (quintic smoothstep). Prevents the throttle
+      ! from amplifying numerical noise in nearly-adiabatic cells.
+
+      ! ::
+
+      superad_reduction_diff_grads_limit = 1d-3
+
+
       ! superad_reduction_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! Implicit alternative to ``okay_to_reduce_gradT_excess``
+      ! Upper saturation cap on the throttle multiplier ``Gamma_factor``.
+      ! The raw ``Gamma_factor`` includes a ``1/sqrt(beta)`` boost that
+      ! can blow up in radiation-dominated cells; this control imposes
+      ! a logistic squash so ``Gamma_factor -> reduction_limit`` as the
+      ! raw value diverges. Set ``-1d0`` to disable the cap (legacy
+      ! behavior); finite values guard against the second MLT call
+      ! seeing ``gradr_scaled`` jumping all the way to ``gradL`` in a
+      ! single Newton iteration.
 
       ! ::
 
-      use_superad_reduction = .false.
-      superad_reduction_Gamma_limit = 0.5d0
-      superad_reduction_Gamma_limit_scale = 5d0
-      superad_reduction_Gamma_inv_scale = 5d0
-      superad_reduction_diff_grads_limit = 1d-3
-      superad_reduction_limit = -1d0
+      superad_reduction_limit = 100d0
+
+
+      ! superad_reduction_use_turnover_limit
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Optional convective-turnover-time limiter for the throttle.
+      ! When ``.true.``, relaxes the throttle using a fraction
+      ! ``f_tau`` set by ``superad_reduction_turnover_limit_function``,
+      ! where ``tau_conv = scale_height(k) / mlt_vc_old(k)`` is the
+      ! local convective turnover time computed from the previous
+      ! step's converged convective velocity. This reflects the
+      ! physical constraint that convection cannot redistribute the
+      ! radiative input on timescales shorter than its own turnover
+      ! time --- so when the global time step ``dt`` is much shorter
+      ! than ``tau_conv`` the throttle should not act at full strength.
+      !
+      ! - ``dt >> tau_conv``: ``f_tau -> 1``, current behavior recovered.
+      ! - ``dt << tau_conv``: ``f_tau -> 0``, throttle smoothly off.
+      !
+      ! ``mlt_vc_old`` is used (rather than the current iterate) so
+      ! that ``tau_conv`` is a frozen parameter w.r.t. the inner
+      ! Newton solve: the auto-diff Jacobian remains clean. When the
+      ! cell was radiative in the previous step
+      ! (``mlt_vc_old = 0``), ``f_tau -> 0`` and the throttle is off
+      ! for this step.
+      !
+      ! Default ``.false.``: the limiter block is skipped and the
+      ! routine behaves exactly as it does without this control.
+
+      ! ::
+
+      superad_reduction_use_turnover_limit = .false.
+
+      ! superad_reduction_turnover_limit_function
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Functional form for the turnover-time relaxation fraction. Allowed
+      ! values are:
+      !
+      ! - ``'exponential'``: ``f_tau = 1 - exp(-dt/tau_conv)``
+      ! - ``'linear'``: ``f_tau = min(dt/tau_conv, 1)``
+
+      ! ::
+
+      superad_reduction_turnover_limit_function = 'exponential'
+
+      ! superad_reduction_turnover_vc_floor_frac
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Fractional floor on mlt_vc_old, expressed as a fraction of the
+      ! local face sound speed csound_face(k), used only when
+      ! superad_reduction_use_turnover_limit = .true.
+      !
+      ! The turnover-time limiter uses
+      !   tau_conv(k) = scale_height(k) / mlt_vc_old(k)
+      ! to gauge whether convection has had time to respond to a step.
+      ! Slow-convection cells (e.g., the surface iron-bump zone) have
+      ! tiny mlt_vc_old, hence huge tau_conv, hence f_tau -> 0 and the
+      ! throttle is suppressed there even when the rest of the star is
+      ! evolving normally. That can leave radiation-pressure inversions
+      ! unsmoothed and make the Newton solve hard.
+      !
+      ! Setting this control to a small positive value floors
+      ! mlt_vc_old at frac * csound_face(k) before computing tau_conv.
+      ! Physical interpretation: even in nominally radiative cells
+      ! perturbations cannot propagate slower than ~ acoustic speed, so
+      ! tau_conv cannot meaningfully exceed scale_height / (frac * cs).
+      ! Numerically this prevents f_tau from dropping arbitrarily close
+      ! to zero in the slow-convection cells most prone to retries.
+      !
+      ! Default 0d0 disables the floor (current behavior). Reasonable
+      ! starting value: 1d-3 (cap tau_conv at ~ scale_height / (1e-3 cs)
+      ! = ~ 1000 sound-crossing times of the scale height).
+
+      ! ::
+
+      superad_reduction_turnover_vc_floor_frac = 0d0
 
 
    ! overshooting

star/defaults/history_columns.list

@@ -935,6 +935,10 @@
       !max_L_rad_div_Ledd
       !max_L_rad_div_Ledd_div_phi_Joss
 
+   !## superad_reduction
+      !max_superad_reduction_factor ! peak Gamma_factor over the star (1 = no reduction)
+      !num_cells_with_superad_reduction ! count of cells where Gamma_factor > 1
+
 
    !## RTI
       !rti_regions <num>

star/defaults/profile_columns.list

@@ -511,6 +511,8 @@
    !dvc_dt_TDC_div_g ! dimensionless ratio of convective velocity to g
 
    !superad_reduction_factor ! gamma_factor from superad_reduction
+   !superad_reduction_Lrad_div_Ledd ! local Lrad/Ledd used by superad_reduction (0 if inactive)
+   !superad_reduction_trigger ! 0=untouched, 1=Eddington, 2=density-inversion, 3=both
    !conv_vel_div_mlt_vc ! conv_vel from any convection model (including RSP) / mlt_vc, for comparison 
 
    !log_Lconv

star/private/adjust_mesh.f90

@@ -410,6 +410,7 @@ integer function remesh(s)
             prv% j_rot, prv% i_rot, prv% omega, prv% D_omega, &
             prv% mlt_vc, prv% lnT, prv% w, specific_PE, specific_KE, &
             prv% m, prv% r, prv% rho, prv% dPdr_dRhodr_info, prv% D_mix, &
+            prv% superad_reduction_factor, &
             cell_type, comes_from, prv% dq, xq_old, s% xh, s% xa, s% dq, xq_new, ierr)
          if (ierr /= 0) then
             s% retry_message = 'do_mesh_adjust failed in mesh_adjust'
@@ -429,6 +430,7 @@ integer function remesh(s)
             s% prev_mesh_omega(k) = prv% prev_mesh_omega(k)
             s% prev_mesh_dq(k) = prv% prev_mesh_dq(k)
             s% prev_mesh_mlt_vc(k) = prv% prev_mesh_mlt_vc(k)
+            s% prev_mesh_superad_reduction_factor(k) = prv% prev_mesh_superad_reduction_factor(k)
          end do
 
          ! restore ST info (for time smoothing)

star/private/alloc.f90

@@ -383,6 +383,8 @@ subroutine star_info_old_arrays(s, action, ierr)
          if (failed('dq_old')) return
          call do1D(s, s% mlt_vc_old, nz, action, ierr)
          if (failed('mlt_vc_old')) return
+         call do1D(s, s% superad_reduction_factor_old, nz, action, ierr)
+         if (failed('superad_reduction_factor_old')) return
          call do1D(s, s% omega_old, nz, action, ierr)
          if (failed('omega_old')) return
          call do1D(s, s% j_rot_old, nz, action, ierr)
@@ -776,6 +778,10 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             if (failed('gradT_excess_effect')) exit
             call do1(s% superad_reduction_factor, c% superad_reduction_factor)
             if (failed('superad_reduction_factor')) exit
+            call do1(s% superad_reduction_Lrad_div_Ledd, c% superad_reduction_Lrad_div_Ledd)
+            if (failed('superad_reduction_Lrad_div_Ledd')) exit
+            call do1_integer(s% superad_reduction_trigger, c% superad_reduction_trigger)
+            if (failed('superad_reduction_trigger')) exit
 
             call do1(s% domega_dlnR, c% domega_dlnR)
             if (failed('domega_dlnR')) exit
@@ -1397,6 +1403,8 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             if (failed('prev_mesh_omega')) exit
             call do1(s% prev_mesh_mlt_vc, c% prev_mesh_mlt_vc)
             if (failed('prev_mesh_mlt_vc')) exit
+            call do1(s% prev_mesh_superad_reduction_factor, c% prev_mesh_superad_reduction_factor)
+            if (failed('prev_mesh_superad_reduction_factor')) exit
             call do1(s% prev_mesh_dq, c% prev_mesh_dq)
             if (failed('prev_mesh_dq')) exit
             ! These are needed for time-smoothing of ST mixing

star/private/ctrls_io.f90

@@ -136,6 +136,8 @@ module ctrls_io
    gradT_excess_min_center_he4, gradT_excess_max_logT, gradT_excess_min_log_tau_full_on, gradT_excess_max_log_tau_full_off, &
     use_superad_reduction, superad_reduction_gamma_limit, superad_reduction_gamma_limit_scale, D_mix_zero_region_top_q, &
     superad_reduction_gamma_inv_scale, superad_reduction_diff_grads_limit, superad_reduction_limit, &
+    superad_reduction_use_turnover_limit, superad_reduction_turnover_limit_function, &
+    superad_reduction_turnover_vc_floor_frac, &
     make_gradr_sticky_in_solver_iters, min_logT_for_make_gradr_sticky_in_solver_iters, &
     max_logT_for_mlt, thermohaline_coeff, thermohaline_option, mixing_length_alpha, remove_small_D_limit, &
     alt_scale_height_flag, Henyey_MLT_y_param, Henyey_MLT_nu_param, no_MLT_below_shock, mlt_make_surface_no_mixing, &
@@ -1078,6 +1080,9 @@ subroutine store_controls(s, ierr)
  s% superad_reduction_gamma_inv_scale = superad_reduction_gamma_inv_scale
  s% superad_reduction_diff_grads_limit = superad_reduction_diff_grads_limit
  s% superad_reduction_limit = superad_reduction_limit
+ s% superad_reduction_use_turnover_limit = superad_reduction_use_turnover_limit
+ s% superad_reduction_turnover_limit_function = superad_reduction_turnover_limit_function
+ s% superad_reduction_turnover_vc_floor_frac = superad_reduction_turnover_vc_floor_frac
 
  s% max_logT_for_mlt = max_logT_for_mlt
  s% mlt_make_surface_no_mixing = mlt_make_surface_no_mixing
@@ -2801,6 +2806,9 @@ subroutine set_controls_for_writing(s, ierr)
  superad_reduction_gamma_inv_scale = s% superad_reduction_gamma_inv_scale
  superad_reduction_diff_grads_limit = s% superad_reduction_diff_grads_limit
  superad_reduction_limit = s% superad_reduction_limit
+ superad_reduction_use_turnover_limit = s% superad_reduction_use_turnover_limit
+ superad_reduction_turnover_limit_function = s% superad_reduction_turnover_limit_function
+ superad_reduction_turnover_vc_floor_frac = s% superad_reduction_turnover_vc_floor_frac
 
  max_logT_for_mlt = s% max_logT_for_mlt
  mlt_make_surface_no_mixing = s% mlt_make_surface_no_mixing

star/private/evolve.f90

@@ -1864,6 +1864,7 @@ integer function prepare_for_new_step(s)
                s% prev_mesh_omega(k) = s% omega(k)
                s% prev_mesh_dq(k) = s% dq(k)
                s% prev_mesh_mlt_vc(k) = s% mlt_vc(k)
+               s% prev_mesh_superad_reduction_factor(k) = s% superad_reduction_factor(k)
                s% prev_mesh_species_or_nvar_hydro_changed = .false.
             end do
             s% prev_mesh_nz = s% nz
@@ -1888,6 +1889,21 @@ integer function prepare_for_new_step(s)
             end if
          end if
 
+         ! Snapshot s%superad_reduction_factor into _old here -- AFTER mesh
+         ! adjust has remapped it onto the new mesh, but BEFORE
+         ! set_vars_if_needed below would overwrite it with the new step's
+         ! first-pass value. Used by the smooth-relaxation turnover-time
+         ! limiter as the Gamma_factor_old anchor.
+         if (s% use_superad_reduction .and. .not. s% RSP_flag) then
+            if (associated(s% superad_reduction_factor) .and. &
+                associated(s% superad_reduction_factor_old)) then
+               do k = 1, s% nz
+                  s% superad_reduction_factor_old(k) = s% superad_reduction_factor(k)
+               end do
+               s% have_superad_reduction_factor = .true.
+            end if
+         end if
+
          call set_vars_if_needed(s, s% dt_next, 'prepare_for_new_step', ierr)
          if (failed('set_vars_if_needed ierr')) return
          call set_phot_info(s)  ! this sets Teff and other stuff
@@ -2198,6 +2214,7 @@ integer function prepare_to_retry(id)
                   s% omega_old(k) = s% prev_mesh_omega(k)
                   s% j_rot_old(k) = s% prev_mesh_j_rot(k)
                   s% mlt_vc_old(k) = s% prev_mesh_mlt_vc(k)
+                  s% superad_reduction_factor_old(k) = s% prev_mesh_superad_reduction_factor(k)
                end do
                call normalize_dqs(s, s% prev_mesh_nz, s% dq_old, ierr)
                if (ierr /= 0) then