superad turnover time limiter

docs/source/changelog.rst

@@ -41,6 +41,8 @@ GYRE has been upgraded to 9.1.1, the most recent stable release. Changes since t
 
 The `MESA SDK <http://user.astro.wisc.edu/~townsend/static.php?ref=mesasdk>`__ recommended for compiling MESA has been updated to 26.6.1. Although this newer SDK is not required to successfully build MESA, it brings the benefit of restoring the cross-platform bit-for-bit compatibility that MESA once enjoyed (meaning that runs on Linux/Intel, MacOS/Intel and MacOS/ARM give identical results).
 
+The new control ``superad_reduction_use_turnover_limit`` relaxes the applied superadiabatic reduction from its previous accepted value toward the instantaneous value. ``superad_reduction_turnover_limit_function`` selects either the exponential response ``1-exp(-dt/tau_conv)`` or the linear response ``min(dt/tau_conv,1)``. The limiter acts on the applied reduction ``1/Gamma_factor``. Zones with a lagged convective velocity use ``scale_height/max(mlt_vc,1d-10 cm/s)``; other zones use the Brunt frequency. The timescale is set at the start of the step and held fixed during solver iterations. With ``use_face_reconstruction``, this calculation uses the reconstructed face thermodynamic state. The scale height is the interpolated or reconstructed face value used by MLT and TDC. The previous reduction is preserved across retries, remeshing, and photo restarts. For ``k > 0``, ``superad_reduction_max_logT`` restricts the reduction to faces whose start-of-step temperature is below the selected ``logT``. Its default is ``7d0``, corresponding to :math:`10^7\,\mathrm{K}`. The ``k=0`` model-construction path is unchanged.
+
 Additional controls are available for TDC envelope remeshing:
 
 - ``remesh_for_TDC_pulsations_when_load`` remeshes a model after loading it from a ``.mod`` file.
@@ -68,6 +70,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`` used ``superad_reduction_Gamma_limit_scale`` instead of ``superad_reduction_Gamma_inv_scale``. This made the latter control inactive. The published algorithm is restored (Jermyn et al. 2023, equation 64, ``alpha_2`` term). Existing test suite cases set the two scales to the same value and are bit-for-bit unaffected. The fix changes results only when the scales differ.
+
+Photo restarts now preserve the stored superadiabatic turnover-limiter state while ``finish_load_model`` rebuilds derived quantities. Previously, the stored reduction could be reset before ``extras_startup`` restored controls selected dynamically by ``run_star_extras``. Runs that leave superadiabatic reduction disabled still reset the factor on the first evolution step.
+
 .. 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``).
 
 Other Changes

linters/check_photos.py

@@ -54,6 +54,9 @@ def discard(self, value):
     "net_iso",
     "model_profile_filename",
     "include_binary_history_in_log_file",
+    "superad_reduction_factor_old",  # Initialized from superad_reduction_factor
+    "have_superad_reduction_factor",  # Set when superad_reduction_factor is read
+    "okay_to_set_superad_reduction_factor",  # Reset during photo loading
     "i_dj_rot_dt",  # Shadows s% i_j_rot
     "i_equ_w_div_wc",  # Shadows s% i_w_div_wc
     "have_j_rot",  # Set from s% rotation_flag

star/defaults/controls.defaults

@@ -2658,29 +2658,113 @@
 
       ! use_superad_reduction
       ! ~~~~~~~~~~~~~~~~~~~~~
+
+      ! Implicit alternative to ``okay_to_reduce_gradT_excess``. MLT or TDC
+      ! is repeated with ``gradr`` moved toward ``gradL`` in cells near the
+      ! Eddington or density-inversion limits. The two contributions are
+      ! evaluated separately and added to ``Gamma_factor``, with
+      ! ``gradr_scaled = gradL + (gradr-gradL)/Gamma_factor``.
+
+      ! ::
+
+      use_superad_reduction = .false.
+
+
+      ! superad_reduction_max_logT
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! For k > 0, apply superadiabatic reduction only where the start-of-step
+      ! ``logT`` is less than this limit. The start-of-step value keeps the
+      ! selection fixed during solver iterations.
+
+      ! ::
+
+      superad_reduction_max_logT = 7d0
+
+
       ! superad_reduction_Gamma_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Threshold in ``Lrad/Ledd`` for the Eddington contribution.
+
+      ! ::
+
+      superad_reduction_Gamma_limit = 0.5d0
+
+
       ! superad_reduction_Gamma_limit_scale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Scale factor for the Eddington contribution. Set to zero to disable
+      ! this contribution.
+
+      ! ::
+
+      superad_reduction_Gamma_limit_scale = 5d0
+
+
       ! superad_reduction_Gamma_inv_scale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Scale factor for the density-inversion contribution, which activates
+      ! when ``Lrad/Ledd > 4*(1-beta)/(4-3*beta)``. Set to zero to disable
+      ! this contribution.
+
+      ! ::
+
+      superad_reduction_Gamma_inv_scale = 5d0
+
+
       ! superad_reduction_diff_grads_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Smoothly turns off the reduction as ``gradT - gradL`` approaches
+      ! this value.
+
+      ! ::
+
+      superad_reduction_diff_grads_limit = 1d-3
+
+
       ! superad_reduction_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! Implicit alternative to ``okay_to_reduce_gradT_excess``
+      ! Upper limit for ``Gamma_factor``. A value less than or equal to one
+      ! disables the limit.
 
       ! ::
 
-      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_use_turnover_limit
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Relaxes the applied reduction from its previous accepted value toward
+      ! the instantaneous value over the local response time. The response
+      ! fraction is set by ``superad_reduction_turnover_limit_function``.
+      ! In zones with a lagged convective velocity, ``tau_conv`` is
+      ! ``scale_height/max(mlt_vc,1d-10 cm/s)``. Other zones use the Brunt
+      ! frequency. The timescale is set at the start of the step and held fixed
+      ! during solver iterations. The limiter acts on ``1/Gamma_factor``, which
+      ! multiplies ``gradr-gradL``, and uses the same interpolated or
+      ! reconstructed face state as MLT or TDC.
+
+      ! ::
+
+      superad_reduction_use_turnover_limit = .false.
+
+      ! superad_reduction_turnover_limit_function
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Sets the response fraction used by the turnover limiter. The
+      ! ``'exponential'`` option uses ``1-exp(-dt/tau_conv)`` and the
+      ! ``'linear'`` option uses ``min(dt/tau_conv,1)``.
+
+      ! ::
+
+      superad_reduction_turnover_limit_function = 'linear'
+
    ! overshooting
    ! ____________
 

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,8 @@ 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/adjust_mesh_split_merge.f90

@@ -657,6 +657,7 @@ subroutine do_merge(s, i_merge, species, new_xa, ierr)
             s% lnT(im) = s% lnT(i0)
             s% D_mix(im) = s% D_mix(i0)
             s% mlt_vc(im) = s% mlt_vc(i0)
+            s% superad_reduction_factor(im) = s% superad_reduction_factor(i0)
             s% csound(im) = s% csound(i0)
             s% tau(im) = s% tau(i0)
             s% opacity(im) = s% opacity(i0)
@@ -831,7 +832,8 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr)
             sumx, sumxp, new_xaL, new_xaR, star_PE0, star_PE1, &
             grad_alpha, f, new_alphaL, new_alphaR, v_R, v_C, v_L, min_dm, &
             mlt_vcL, mlt_vcR, tauL, tauR, etrb, etrb_L, etrb_C, etrb_R, grad_etrb, &
-            j_rot_new, dmbar_old, dmbar_p1_old, dmbar_new, dmbar_p1_new, dmbar_p2_new, J_old
+            j_rot_new, dmbar_old, dmbar_p1_old, dmbar_new, dmbar_p1_new, dmbar_p2_new, J_old, &
+            superad_reduction_factorL, superad_reduction_factorR
          logical :: done, use_new_grad_rho
          include 'formats'
 
@@ -868,13 +870,16 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr)
 
          rR = s% r(i)
          mlt_vcR = s% mlt_vc(i)
+         superad_reduction_factorR = s% superad_reduction_factor(i)
          if (i == nz) then
             rL = s% R_center
             mlt_vcL = 0d0
+            superad_reduction_factorL = 1d0
             tauL = tau_center
          else
             rL = s% r(ip)
             mlt_vcL = s% mlt_vc(ip)
+            superad_reduction_factorL = s% superad_reduction_factor(ip)
             tauL = s% tau(ip)
          end if
 
@@ -1038,6 +1043,7 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr)
                s% lnT(jp) = s% lnT(j)
                s% D_mix(jp) = s% D_mix(j)
                s% mlt_vc(jp) = s% mlt_vc(j)
+               s% superad_reduction_factor(jp) = s% superad_reduction_factor(j)
                s% csound(jp) = s% csound(j)
                s% tau(jp) = s% tau(j)
                s% opacity(jp) = s% opacity(j)
@@ -1175,9 +1181,11 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr)
             s% dPdr_dRhodr_info(ip) = s% dPdr_dRhodr_info(i)
          end if
 
-         ! mlt_vc is face-based, so a split creates a new interior face value here.
+         ! These are face-based, so a split creates new interior face values here.
          ! Interpolate using the same left/right orientation as tau(ip).
          s% mlt_vc(ip) = mlt_vcR + (mlt_vcL - mlt_vcR)*dMR/dM
+         s% superad_reduction_factor(ip) = superad_reduction_factorR + &
+            (superad_reduction_factorL - superad_reduction_factorR)*dMR/dM
 
          s% tau(ip) = tauR + (tauL - tauR)*dMR/dM
          if (is_bad(s% tau(ip))) then

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,7 @@ 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
+            if (action == do_allocate) s% superad_reduction_factor = 1d0
 
             call do1(s% domega_dlnR, c% domega_dlnR)
             if (failed('domega_dlnR')) exit
@@ -1163,6 +1166,9 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             if (failed('unsmoothed_brunt_B')) exit
             call do1(s% smoothed_brunt_B, c% smoothed_brunt_B)
             if (failed('smoothed_brunt_B')) exit
+            call do1(s% tau_conv_start, c% tau_conv_start)
+            if (failed('tau_conv_start')) exit
+            if (action == do_allocate) s% tau_conv_start = 0d0
 
             call do1(s% RTI_du_diffusion_kick, c% RTI_du_diffusion_kick)
             if (failed('RTI_du_diffusion_kick')) exit
@@ -1427,6 +1433,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

@@ -138,8 +138,10 @@ module ctrls_io
    gradT_excess_beta1, gradT_excess_lambda2, gradT_excess_beta2, &
    gradT_excess_dlambda, gradT_excess_dbeta, gradT_excess_max_center_h1, &
    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, &
+    use_superad_reduction, superad_reduction_max_logT, &
+    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, &
     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, &
@@ -656,6 +658,16 @@ subroutine check_controls(s, ierr)
        return
     end if
 
+    if (.not. (trim(s% superad_reduction_turnover_limit_function) == 'linear' .or. &
+          trim(s% superad_reduction_turnover_limit_function) == 'exponential')) then
+       write(*,'(A)')
+       write(*,*) "Invalid choice for superad_reduction_turnover_limit_function"
+       write(*,*) "Available options are 'linear' or 'exponential'"
+       write(*,'(A)')
+       ierr = -1
+       return
+    end if
+
  end subroutine check_controls
 
 
@@ -1041,11 +1053,14 @@ subroutine store_controls(s)
  s% dq_D_mix_zero_at_H_C_crossover = dq_D_mix_zero_at_H_C_crossover
 
  s% use_superad_reduction = use_superad_reduction
+ s% superad_reduction_max_logT = superad_reduction_max_logT
  s% superad_reduction_gamma_limit = superad_reduction_gamma_limit
  s% superad_reduction_gamma_limit_scale = superad_reduction_gamma_limit_scale
  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% max_logT_for_mlt = max_logT_for_mlt
  s% mlt_make_surface_no_mixing = mlt_make_surface_no_mixing
@@ -2769,11 +2784,14 @@ subroutine set_controls_for_writing(s, ierr)
  dq_D_mix_zero_at_H_C_crossover = s% dq_D_mix_zero_at_H_C_crossover
 
  use_superad_reduction = s% use_superad_reduction
+ superad_reduction_max_logT = s% superad_reduction_max_logT
  superad_reduction_gamma_limit = s% superad_reduction_gamma_limit
  superad_reduction_gamma_limit_scale = s% superad_reduction_gamma_limit_scale
  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
 
  max_logT_for_mlt = s% max_logT_for_mlt
  mlt_make_surface_no_mixing = s% mlt_make_surface_no_mixing

star/private/evolve.f90

@@ -319,6 +319,7 @@ integer function do_step_part1(id, first_try)
          s% need_to_setvars = .true.  ! always start fresh
          s% okay_to_set_mixing_info = .true.  ! set false by element diffusion
          s% okay_to_set_mlt_vc = .false.  ! don't change mlt_vc until have set mlt_vc_old
+         s% okay_to_set_superad_reduction_factor = .false.
 
          if (s% timestep_hold > s% model_number + 10000) then
             write(*,3) 'ERROR: s% timestep_hold', s% timestep_hold, s% model_number
@@ -698,6 +699,7 @@ integer function do_step_part2(id, first_try)
                s% have_mlt_vc = .true.
             end if
             s% okay_to_set_mlt_vc = .false.
+            s% okay_to_set_superad_reduction_factor = .false.
          end if
 
          if (.not. okay_energy_conservation()) return
@@ -1864,6 +1866,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
@@ -2198,6 +2201,8 @@ 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

star/private/evolve_support.f90

@@ -58,6 +58,12 @@ subroutine new_generation(s, ierr)
             call copy_to_old(s% mlt_vc, s% mlt_vc_old, ierr)
             if (ierr /= 0) return
 
+            call copy_to_old(s% superad_reduction_factor, &
+               s% superad_reduction_factor_old, ierr)
+            if (ierr /= 0) return
+            s% have_superad_reduction_factor = .true.
+            s% okay_to_set_superad_reduction_factor = .true.
+
             call enlarge_if_needed_2(s% xh_old,s% nvar_hydro,nz,nz_alloc_extra,ierr)
             if (ierr /= 0) return
             if (s% fill_arrays_with_NaNs) call fill_with_NaNs_2d(s% xh_old)
@@ -191,8 +197,10 @@ subroutine set_current_to_old(s)
                end do
                s% dq(k) = s% dq_old(k)
                s% mlt_vc(k) = s% mlt_vc_old(k)
+               s% superad_reduction_factor(k) = s% superad_reduction_factor_old(k)
             end do
             s% okay_to_set_mlt_vc = .true.
+            s% okay_to_set_superad_reduction_factor = .true.
 
             call set_qs(s, s% nz, s% q, s% dq, ierr)
             if (ierr /= 0) then