From 15a32f5c4759e0f27708b569ecc8f3c041b7ccec Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Tue, 19 May 2026 10:36:14 -0400 Subject: [PATCH 01/14] star/turb_support: fix dead-code in superad_reduction inversion contribution In set_superad_reduction, the second contribution to Gamma_term (the Joss--Salpeter--Ostriker density-inversion criterion, Paxton+2013 eq. 17) was scaled by superad_reduction_Gamma_limit_scale instead of the intended superad_reduction_Gamma_inv_scale. This made the Gamma_inv_scale control silently inert: changing it had no effect on the inversion-driven throttle. Two-line fix: scale_value1 -> scale_value2 in both branches of the inversion contribution. Restores the algorithm to its as-published form (Jermyn et al. 2023, eq. 64, terms alpha_1 vs alpha_2). Existing test_suite cases all set Gamma_limit_scale = Gamma_inv_scale (both = 5d0, the default), so this change is bit-for-bit identical on every shipped test. The fix is detectable only when the two scales are set to different values. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/private/turb_support.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 67eda0362..f820ea59f 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -494,9 +494,9 @@ subroutine set_superad_reduction() if (Lrad_div_Ledd% val > Gamma_inv_threshold) then alfa0 = Lrad_div_Ledd/Gamma_inv_threshold-1d0 if (alfa0 < 1d0) then - Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0) + Gamma_term = Gamma_term + scale_value2*(0.5d0*alfa0*alfa0) else - Gamma_term = Gamma_term + scale_value1*(alfa0-0.5d0) + Gamma_term = Gamma_term + scale_value2*(alfa0-0.5d0) end if !Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0) end if From f93f5fd385028f7ef9aa6a0d6e400ba4a62a6b63 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Tue, 19 May 2026 10:37:15 -0400 Subject: [PATCH 02/14] star/controls.defaults: cap superad_reduction_limit at 100 by default Two changes: 1. Default for superad_reduction_limit 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, stressing the outer structure solver. The cap acts as a numerical guard rail without changing the underlying physics. Existing test_suite cases that exercise superad_reduction all explicitly set this control to -1d0, so they are unaffected by the default change. 2. Filled in per-control docstrings for the six existing superad_reduction_* controls (previously they all shared a single stub block). Each control now has its own ~~~~~ section explaining its physical role -- referencing Jermyn et al. 2023 eq. 64 and Paxton, Cantiello et al. 2013 eq. 17 (Joss+1973 density-inversion criterion). No code change beyond the cap default. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/defaults/controls.defaults | 78 ++++++++++++++++++++++++++++++--- 1 file changed, 71 insertions(+), 7 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 357fe8631..cae38f1e9 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2655,27 +2655,91 @@ ! 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 ! overshooting From 040250ab7adc5f5b01f76fe95a1ff663b53d890e Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Tue, 19 May 2026 10:41:16 -0400 Subject: [PATCH 03/14] star/superad_reduction: per-cell + history diagnostics Adds three new diagnostic outputs for use_superad_reduction so users can see what the throttle is doing rather than just observing its downstream effects on gradT. Profile columns (opt-in, in star/defaults/profile_columns.list): superad_reduction_factor (already existed, kept) superad_reduction_Lrad_div_Ledd ! L_rad/L_Edd computed inside ! set_superad_reduction from the ! first-pass (unthrottled) gradT superad_reduction_trigger ! 0 = throttle inactive ! 1 = activated by Eddington proximity ! (Lrad/Ledd > Gamma_limit) ! 2 = activated by density-inversion ! criterion alone ! (Joss+1973 = Paxton+2013 eq. 17) ! 3 = both triggers active History columns (opt-in, in star/defaults/history_columns.list): max_superad_reduction_factor (already existed, kept) num_cells_with_superad_reduction ! count of cells where factor > 1 Wired through alloc.f90, profile_getval.f90, star_profile_def.f90, history.f90, and star_history_def.f90; new pointer arrays declared in star_data/public/star_data_step_work.inc. No physics change -- the throttle algorithm is untouched, only its per-cell state is exposed. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/defaults/history_columns.list | 4 ++++ star/defaults/profile_columns.list | 2 ++ star/private/alloc.f90 | 4 ++++ star/private/history.f90 | 16 ++++++++++++++++ star/private/profile_getval.f90 | 5 +++++ star/private/star_history_def.f90 | 8 +++++++- star/private/star_profile_def.f90 | 6 +++++- star/private/turb_support.f90 | 13 ++++++++++++- star_data/public/star_data_step_work.inc | 13 +++++++++++++ 9 files changed, 68 insertions(+), 3 deletions(-) diff --git a/star/defaults/history_columns.list b/star/defaults/history_columns.list index 858f88e3d..06fad8c28 100644 --- a/star/defaults/history_columns.list +++ b/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 diff --git a/star/defaults/profile_columns.list b/star/defaults/profile_columns.list index 2bf7338c2..4b683c4c8 100644 --- a/star/defaults/profile_columns.list +++ b/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 diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 675345b21..3de795202 100644 --- a/star/private/alloc.f90 +++ b/star/private/alloc.f90 @@ -776,6 +776,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 diff --git a/star/private/history.f90 b/star/private/history.f90 index 66540268e..3731311e7 100644 --- a/star/private/history.f90 +++ b/star/private/history.f90 @@ -1993,6 +1993,22 @@ subroutine history_getval(& case(h_gradT_excess_max_lambda) val = s% gradT_excess_max_lambda + case(h_max_superad_reduction_factor) + val = 1d0 + if (s% use_superad_reduction) then + do k = 1, nz + if (s% superad_reduction_factor(k) > val) val = s% superad_reduction_factor(k) + end do + end if + case(h_num_cells_with_superad_reduction) + int_val = 0 + if (s% use_superad_reduction) then + do k = 1, nz + if (s% superad_reduction_factor(k) > 1d0) int_val = int_val + 1 + end do + end if + is_int_val = .true. + case(h_max_L_rad_div_Ledd) do k = 1, nz tmp = get_Lrad_div_Ledd(s, k) diff --git a/star/private/profile_getval.f90 b/star/private/profile_getval.f90 index 471399f94..9e2cf4bfe 100644 --- a/star/private/profile_getval.f90 +++ b/star/private/profile_getval.f90 @@ -499,6 +499,11 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) case (p_superad_reduction_factor) val = s% superad_reduction_factor(k) + case (p_superad_reduction_Lrad_div_Ledd) + val = s% superad_reduction_Lrad_div_Ledd(k) + case (p_superad_reduction_trigger) + int_val = s% superad_reduction_trigger(k) + int_flag = .true. case (p_gradT_excess_effect) val = s% gradT_excess_effect(k) case (p_diff_grads) diff --git a/star/private/star_history_def.f90 b/star/private/star_history_def.f90 index e569989d9..b87f2426d 100644 --- a/star/private/star_history_def.f90 +++ b/star/private/star_history_def.f90 @@ -402,7 +402,10 @@ module star_history_def integer, parameter :: h_gradT_excess_max_lambda = h_gradT_excess_min_beta + 1 integer, parameter :: h_gradT_excess_alpha = h_gradT_excess_max_lambda + 1 - integer, parameter :: h_log_Ledd = h_gradT_excess_alpha + 1 + integer, parameter :: h_max_superad_reduction_factor = h_gradT_excess_alpha + 1 + integer, parameter :: h_num_cells_with_superad_reduction = h_max_superad_reduction_factor + 1 + + integer, parameter :: h_log_Ledd = h_num_cells_with_superad_reduction + 1 integer, parameter :: h_compactness = h_log_Ledd + 1 integer, parameter :: h_compactness_parameter = h_compactness + 1 integer, parameter :: h_mu4 = h_compactness_parameter + 1 @@ -915,6 +918,9 @@ subroutine history_column_names_init(ierr) history_column_name(h_gradT_excess_min_beta) = 'gradT_excess_min_beta' history_column_name(h_gradT_excess_max_lambda) = 'gradT_excess_max_lambda' + history_column_name(h_max_superad_reduction_factor) = 'max_superad_reduction_factor' + history_column_name(h_num_cells_with_superad_reduction) = 'num_cells_with_superad_reduction' + history_column_name(h_gamma1_min) = 'gamma1_min' history_column_name(h_logT_max) = 'logT_max' history_column_name(h_logQ_max) = 'logQ_max' diff --git a/star/private/star_profile_def.f90 b/star/private/star_profile_def.f90 index a51dd25c6..bfd722683 100644 --- a/star/private/star_profile_def.f90 +++ b/star/private/star_profile_def.f90 @@ -72,7 +72,9 @@ module star_profile_def integer, parameter :: p_vel_km_per_s = p_v_kms + 1 integer, parameter :: p_log_abs_v = p_vel_km_per_s + 1 integer, parameter :: p_superad_reduction_factor = p_log_abs_v + 1 - integer, parameter :: p_gradT_excess_effect = p_superad_reduction_factor + 1 + integer, parameter :: p_superad_reduction_Lrad_div_Ledd = p_superad_reduction_factor + 1 + integer, parameter :: p_superad_reduction_trigger = p_superad_reduction_Lrad_div_Ledd + 1 + integer, parameter :: p_gradT_excess_effect = p_superad_reduction_trigger + 1 integer, parameter :: p_log_diff_grads = p_gradT_excess_effect + 1 integer, parameter :: p_diff_grads = p_log_diff_grads + 1 @@ -773,6 +775,8 @@ subroutine profile_column_names_init(ierr) profile_column_name(p_diff_grads) = 'diff_grads' profile_column_name(p_gradT_excess_effect) = 'gradT_excess_effect' profile_column_name(p_superad_reduction_factor) = 'superad_reduction_factor' + profile_column_name(p_superad_reduction_Lrad_div_Ledd) = 'superad_reduction_Lrad_div_Ledd' + profile_column_name(p_superad_reduction_trigger) = 'superad_reduction_trigger' profile_column_name(p_v) = 'v' profile_column_name(p_velocity) = 'velocity' diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index f820ea59f..7ab102e7e 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -287,7 +287,11 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg conv_vel = 0d0 D = 0d0 Gamma = 0d0 - if (k /= 0) s% superad_reduction_factor(k) = 1d0 + if (k /= 0) then + s% superad_reduction_factor(k) = 1d0 + s% superad_reduction_Lrad_div_Ledd(k) = 0d0 + s% superad_reduction_trigger(k) = 0 + end if ! Bail if we asked for no mixing, or if parameters are bad. if (MLT_option == 'none' .or. beta < 1d-10 .or. mixing_length_alpha <= 0d0 .or. & @@ -465,6 +469,11 @@ subroutine set_superad_reduction() Lrad_div_Ledd = 4d0*crad/3d0*pow4(T)/P*gradT Gamma_inv_threshold = 4d0*(1d0-beta)/(4d0-3*beta) + if (k /= 0) then + s% superad_reduction_Lrad_div_Ledd(k) = Lrad_div_Ledd% val + s% superad_reduction_trigger(k) = 0 + end if + Gamma_factor = 1d0 if (gradT > gradL) then if (Lrad_div_Ledd > Gamma_limit .or. Lrad_div_Ledd > Gamma_inv_threshold) then @@ -483,6 +492,7 @@ subroutine set_superad_reduction() ! Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0) !end if if (Lrad_div_Ledd > Gamma_limit) then + if (k /= 0) s% superad_reduction_trigger(k) = s% superad_reduction_trigger(k) + 1 alfa0 = Lrad_div_Ledd/Gamma_limit-1d0 if (alfa0 < 1d0) then Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0) @@ -492,6 +502,7 @@ subroutine set_superad_reduction() !Gamma_term = Gamma_term + scale_value1*pow2(Lrad_div_Ledd/Gamma_limit-1d0) end if if (Lrad_div_Ledd% val > Gamma_inv_threshold) then + if (k /= 0) s% superad_reduction_trigger(k) = s% superad_reduction_trigger(k) + 2 alfa0 = Lrad_div_Ledd/Gamma_inv_threshold-1d0 if (alfa0 < 1d0) then Gamma_term = Gamma_term + scale_value2*(0.5d0*alfa0*alfa0) diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 1991b965d..c69242038 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -355,6 +355,19 @@ real(dp), pointer :: superad_reduction_factor(:) ! only set when use_superad_reduction = .true. ! reports effect of superad_reduction. value = 1d0 for no effect. + real(dp), pointer :: superad_reduction_Lrad_div_Ledd(:) + ! only set when use_superad_reduction = .true. + ! reports Lrad/Ledd = (4 a T^4 / 3 P) * gradT, computed from the + ! first-pass (unthrottled) gradT inside set_superad_reduction. + ! Diagnostic only; not used by the solver. + + integer, pointer :: superad_reduction_trigger(:) + ! only set when use_superad_reduction = .true. + ! 0 = throttle inactive in this cell. + ! 1 = activated by Eddington threshold (Lrad/Ledd > Gamma_limit). + ! 2 = activated by density-inversion threshold only. + ! 3 = both triggers active. + ! mixing length alpha for MLT (can be set by user and can vary within model) real(dp), pointer :: alpha_mlt(:) From 248e1d4b44f4ed269bd9bec86fbca1f0f11de8d0 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Tue, 19 May 2026 10:52:56 -0400 Subject: [PATCH 04/14] docs/changelog: document superad_reduction A1 fix, A2 cap default, A3 diagnostics Three bullets under "Changes in main": - New Features: the new superad_reduction_Lrad_div_Ledd / superad_reduction_trigger profile columns and the num_cells_with_superad_reduction history column. - Bug Fixes (A1): scale_value1 -> scale_value2 in the inversion contribution to Gamma_term. Restores the algorithm to its as-published form per Jermyn et al. 2023 eq. 64. Bit-for-bit unchanged on existing test_suite cases. - Bug Fixes (A2): default superad_reduction_limit changed -1d0 -> 100d0. Activates the documented cap that was previously off by default. Bit-for-bit unchanged on existing test_suite cases. Co-Authored-By: Claude Opus 4.7 (1M context) --- docs/source/changelog.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 74fc7d8a3..498a50af6 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -30,6 +30,8 @@ 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. + .. _Bug Fixes main: Bug Fixes @@ -45,6 +47,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``). From 7c42bf508b3a0ffcf858db141cf9868f012017d6 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Wed, 20 May 2026 15:40:05 -0400 Subject: [PATCH 05/14] star/superad_reduction: opt-in convective-turnover-time limiter When the new control superad_reduction_use_turnover_limit is .true., multiply the throttle excess (Gamma_factor - 1) by f_tau = 1 - exp(-dt / tau_conv), tau_conv = scale_height / mlt_vc_old in each convective cell. This is the natural first-order response fraction of a relaxing system with characteristic time tau_conv: the fraction of the steady-state correction realized in time dt. The limiter recovers the current behavior in the dt >> tau_conv limit and smoothly suppresses the throttle in the dt << tau_conv limit -- where the implicit-step approximation is asking convection to respond faster than it can. Implementation: - new logical control with default .false. (opt-in) - applied AFTER the existing logistic cap superad_reduction_limit, so the cap remains an upper bound on Gamma_factor and the limiter just linearly interpolates between 1 and the capped value - uses mlt_vc_old (real(dp), no autoDiff partials) so tau_conv is a frozen exogenous parameter w.r.t. the inner Newton solve -- avoiding the bistable feedback loop that broke our inner-Picard experiments - scale_height is auto_diff in the calling context; its partials propagate cleanly through tau_conv and f_turnover so the outer Newton solver sees an analytically-correct linearization Files: - star_data/private/star_controls.inc: new logical field - star/private/ctrls_io.f90: namelist entry + read/write - star/private/turb_support.f90: limiter block + locals - star/defaults/controls.defaults: default + docstring - docs/source/changelog.rst: new-feature entry Default disabled, so every existing test_suite case is bit-for-bit unchanged. Design document: report/turnover_time_limiter.{tex,pdf}. Co-Authored-By: Claude Opus 4.7 (1M context) --- docs/source/changelog.rst | 2 ++ star/defaults/controls.defaults | 35 +++++++++++++++++++++++++++++ star/private/ctrls_io.f90 | 3 +++ star/private/turb_support.f90 | 16 +++++++++++++ star_data/private/star_controls.inc | 1 + 5 files changed, 57 insertions(+) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 498a50af6..68f1b23ad 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -32,6 +32,8 @@ MESA no longer stops when reactions for which special rates are set are not in t 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. + .. _Bug Fixes main: Bug Fixes diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index cae38f1e9..3a33bf1ef 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2742,6 +2742,41 @@ superad_reduction_limit = 100d0 + ! superad_reduction_use_turnover_limit + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Optional convective-turnover-time limiter for the throttle. + ! When ``.true.``, multiplies the throttle excess + ! ``(Gamma_factor - 1)`` by + ! + ! ``f_tau = 1 - exp(-dt / tau_conv)`` + ! + ! 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. + + ! overshooting ! ____________ diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 3973e3970..04ee1efd0 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -136,6 +136,7 @@ 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, & 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 +1079,7 @@ 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% max_logT_for_mlt = max_logT_for_mlt s% mlt_make_surface_no_mixing = mlt_make_surface_no_mixing @@ -2801,6 +2803,7 @@ 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 max_logT_for_mlt = s% max_logT_for_mlt mlt_make_surface_no_mixing = s% mlt_make_surface_no_mixing diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 7ab102e7e..cd819ccb7 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -221,6 +221,8 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg ! these are used by use_superad_reduction real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit + real(dp) :: vc_old_local + type(auto_diff_real_star_order1) :: tau_conv, f_turnover type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, & diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq, mlt_Pturb, Ptot logical :: test_partials, using_TDC, have_Y_face_guess @@ -520,6 +522,20 @@ subroutine set_superad_reduction() exp_limit = exp(-lambda_limit*(Gamma_factor-1d0)) Gamma_factor = 2d0*(reduction_limit-1d0)*(1d0/(1d0+exp_limit)-0.5d0)+1d0 end if + ! Convective-turnover-time limiter: smoothly suppress the throttle + ! when dt < tau_conv = scale_height / mlt_vc_old. Uses the + ! previous-step converged vc (real(dp), no autoDiff partials) + ! and the step's dt, both frozen as parameters w.r.t. the + ! inner Newton solve. f_turnover -> 1 as dt/tau_conv -> infty + ! (current behavior); f_turnover -> 0 as dt/tau_conv -> 0 + ! (throttle off). Opt-in via superad_reduction_use_turnover_limit. + if (s% superad_reduction_use_turnover_limit .and. & + Gamma_factor > 1d0 .and. k > 0) then + vc_old_local = max(s% mlt_vc_old(k), 1d-30) + tau_conv = scale_height / vc_old_local + f_turnover = 1d0 - exp(-s% dt / tau_conv) + Gamma_factor = 1d0 + f_turnover * (Gamma_factor - 1d0) + end if end if end if end if diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index 2d0c29d8d..a930a9681 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -306,6 +306,7 @@ superad_reduction_Gamma_inv_scale, & superad_reduction_diff_grads_limit, & superad_reduction_limit + logical :: superad_reduction_use_turnover_limit ! mixing parameters From 20f2340ed6374bf55f513d9901cb6314361d20de Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Wed, 20 May 2026 16:44:28 -0400 Subject: [PATCH 06/14] star/turb_support: guard turnover limiter against unassociated mlt_vc_old set_superad_reduction is called during the relax / model-construction phase, before s% mlt_vc_old has been allocated and before the first dt is set. The new turnover-time limiter dereferenced s% mlt_vc_old(k) unconditionally, which segfaulted whenever superad_reduction_use_turnover_limit was .true. Gate the limiter on s% have_mlt_vc, associated(s% mlt_vc_old) and s% dt > 0d0 so the block only runs once a previous-step convective velocity exists. Pre-evolution callers fall through with Gamma_factor unchanged (i.e. the unmodified Jermyn23 reduction is applied). Co-Authored-By: Claude Opus 4.7 (1M context) --- star/private/turb_support.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index cd819ccb7..64ea77f43 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -529,8 +529,12 @@ subroutine set_superad_reduction() ! inner Newton solve. f_turnover -> 1 as dt/tau_conv -> infty ! (current behavior); f_turnover -> 0 as dt/tau_conv -> 0 ! (throttle off). Opt-in via superad_reduction_use_turnover_limit. + ! Skip until previous-step mlt_vc has been populated, otherwise + ! s% mlt_vc_old is unassociated and dereferencing it segfaults. if (s% superad_reduction_use_turnover_limit .and. & - Gamma_factor > 1d0 .and. k > 0) then + Gamma_factor > 1d0 .and. k > 0 .and. & + s% have_mlt_vc .and. associated(s% mlt_vc_old) .and. & + s% dt > 0d0) then vc_old_local = max(s% mlt_vc_old(k), 1d-30) tau_conv = scale_height / vc_old_local f_turnover = 1d0 - exp(-s% dt / tau_conv) From 96ee68e6c042431ee10102467d8b002086d73829 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Wed, 20 May 2026 18:35:24 -0400 Subject: [PATCH 07/14] star/superad_reduction: add v_c^old floor for the turnover-time limiter MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New real control `superad_reduction_turnover_vc_floor_frac` floors mlt_vc_old at frac * csound_face(k) before computing tau_conv = scale_height/mlt_vc_old. Slow-convection iron-bump cells otherwise have tiny mlt_vc_old, hence huge tau_conv, hence f_turnover -> 0 and the throttle is fully suppressed there even when the rest of the star is evolving normally. That leaves radiation-pressure inversions unsmoothed and made the Newton solve grind through hundreds of retries. Validated on the 60M envelope-issues benchmark: with frac=1d-3 the model restarted from photos/x00001000 terminates cleanly at xa_central_lower_limit in 1910 steps with 15 retries — matching the limiter-off baseline (1940 steps, 14 retries) and avoiding the 780-retry / 6000-step grind seen with the unfloored limiter on. Default frac=0d0 disables the floor (current behavior preserved bit-for-bit). Co-Authored-By: Claude Opus 4.7 (1M context) --- star/defaults/controls.defaults | 32 +++++++++++++++++++++++++++++ star/private/ctrls_io.f90 | 4 +++- star/private/turb_support.f90 | 10 +++++++-- star_data/private/star_controls.inc | 1 + 4 files changed, 44 insertions(+), 3 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 3a33bf1ef..c5848219c 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2776,6 +2776,38 @@ superad_reduction_use_turnover_limit = .false. + ! 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 ! ____________ diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 04ee1efd0..e12508ed1 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -136,7 +136,7 @@ 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_use_turnover_limit, 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, & @@ -1080,6 +1080,7 @@ subroutine store_controls(s, ierr) 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_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 @@ -2804,6 +2805,7 @@ subroutine set_controls_for_writing(s, ierr) 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_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 diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 64ea77f43..12ae9125a 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -221,7 +221,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg ! these are used by use_superad_reduction real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit - real(dp) :: vc_old_local + real(dp) :: vc_old_local, vc_old_floor type(auto_diff_real_star_order1) :: tau_conv, f_turnover type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, & diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq, mlt_Pturb, Ptot @@ -535,7 +535,13 @@ subroutine set_superad_reduction() Gamma_factor > 1d0 .and. k > 0 .and. & s% have_mlt_vc .and. associated(s% mlt_vc_old) .and. & s% dt > 0d0) then - vc_old_local = max(s% mlt_vc_old(k), 1d-30) + ! Optional floor on mlt_vc_old at a fraction of the local + ! face sound speed. Prevents tau_conv from blowing up in + ! slow-convection iron-bump cells where f_turnover would + ! otherwise drop to ~0 and leave the throttle off. + vc_old_floor = s% superad_reduction_turnover_vc_floor_frac & + * s% csound_face(k) + vc_old_local = max(s% mlt_vc_old(k), vc_old_floor, 1d-30) tau_conv = scale_height / vc_old_local f_turnover = 1d0 - exp(-s% dt / tau_conv) Gamma_factor = 1d0 + f_turnover * (Gamma_factor - 1d0) diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index a930a9681..bd86209b1 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -307,6 +307,7 @@ superad_reduction_diff_grads_limit, & superad_reduction_limit logical :: superad_reduction_use_turnover_limit + real(dp) :: superad_reduction_turnover_vc_floor_frac ! mixing parameters From d2aa8eb4505f87af6cb72fbc5e58947ce7478dd4 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Wed, 20 May 2026 23:06:17 -0400 Subject: [PATCH 08/14] docs/changelog: document v_c^old floor for the turnover limiter Adds a paragraph noting that superad_reduction_turnover_vc_floor_frac exists, what it does (floor mlt_vc_old at frac * csound_face), and why it matters (without the floor, slow-convection iron-bump cells zero out the throttle and the Newton solver bogs down). Records the 60 M_sun validation: 15 retries with floor=1d-3 vs 780 retries without the floor, matching the limiter-off baseline (14 retries). Co-Authored-By: Claude Opus 4.7 (1M context) --- docs/source/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 68f1b23ad..3c079280e 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -34,6 +34,8 @@ Two new profile columns, ``superad_reduction_Lrad_div_Ledd`` and ``superad_reduc 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. + .. _Bug Fixes main: Bug Fixes From 1cece83a0d2f7f428fe02a472496b903ed261508 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Thu, 21 May 2026 19:41:07 -0400 Subject: [PATCH 09/14] star/superad_reduction: relax Gamma_fac smoothly from the previous-step value The previous turnover-limiter formula Gamma_fac_new = 1 + f_tau * (Gamma_fac_inst - 1) relaxes Gamma_fac toward 1 (no throttle) when dt << tau_conv. That is physically wrong: it implies convection can instantaneously release the throttle. The model state in the previous step was adapted to some throttle value Gamma_fac_old; abruptly resetting Gamma_fac to 1 because the new dt is short produces a discontinuous structural rearrangement. In practice this manifests as a small but visible HRD jump at the end of the to_cc phase, when dt collapses by orders of magnitude during iron-core formation. Correct first-order relaxation: anchor at last-step Gamma_fac_old and relax toward the instantaneous value with characteristic time tau_conv, Gamma_fac_new = Gamma_fac_old + f_tau * (Gamma_fac_inst - Gamma_fac_old) = (1 - f_tau) * Gamma_fac_old + f_tau * Gamma_fac_inst When dt >> tau_conv (steady evolution): f_tau -> 1, Gamma_fac -> Gamma_fac_inst. Identical to the previous behaviour. When dt << tau_conv: f_tau -> 0, Gamma_fac -> Gamma_fac_old. Convection had no time to adapt, so the throttle is held at the previous-step value. This is the discrete form of dGamma_fac/dt = (Gamma_fac_inst - Gamma_fac)/tau_conv, the natural relaxation equation for a quantity with response time tau_conv. Plumbing -------- - star_data/public/star_data_step_work.inc: declare superad_reduction_factor_old(:) and have_superad_reduction_factor alongside the other paired _old arrays. - star/private/alloc.f90: allocate superad_reduction_factor_old in star_info_old_arrays so check_sizes does not trip on a null pointer. - star/private/evolve_support.f90: copy_to_old(superad_reduction_factor, superad_reduction_factor_old) at the end of each successful step. Done unconditionally so check_sizes is happy; have_superad_reduction_factor gates whether the limiter actually consumes the snapshot. - star/private/{init,remove_shells}.f90: initialise the flag to false at startup and on model reload. - star/private/turb_support.f90: the formula change. Guard adds s%have_superad_reduction_factor and associated(...) to the existing 5-way guard. Clamp Gamma_fac_old at >= 1 to avoid relaxing from a non-physical anchor. Default-on: the new formula reproduces the old one exactly on the first step (when have_superad_reduction_factor is false) and is strictly more correct on every subsequent step. No new control needed. Validation: 60 M_sun benchmark, to_cc phase, restart from after_core_c_burn. Both old and new formulae reach fe_core_infall_limit and end with identical central T_c / rho_c (Si-burning corner). The old formula's HRD endpoint sits at log Teff=3.58 / log L=5.58 -- the visible jump. The new formula's endpoint sits at log Teff=3.70 / log L=6.11, on a smooth continuation of the track. See report/smooth_vs_old_limiter.pdf. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/private/alloc.f90 | 2 + star/private/evolve_support.f90 | 14 +++++++ star/private/init.f90 | 1 + star/private/remove_shells.f90 | 1 + star/private/turb_support.f90 | 53 +++++++++++++++++------- star_data/public/star_data_step_work.inc | 15 ++++++- 6 files changed, 70 insertions(+), 16 deletions(-) diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 3de795202..33af7dbef 100644 --- a/star/private/alloc.f90 +++ b/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) diff --git a/star/private/evolve_support.f90 b/star/private/evolve_support.f90 index 46c692e92..f463a22d2 100644 --- a/star/private/evolve_support.f90 +++ b/star/private/evolve_support.f90 @@ -58,6 +58,20 @@ subroutine new_generation(s, ierr) call copy_to_old(s% mlt_vc, s% mlt_vc_old, ierr) if (ierr /= 0) return + ! Smooth-relaxation limiter for use_superad_reduction needs + ! the previous step's converged Gamma_factor. Mirror the + ! mlt_vc -> mlt_vc_old pattern. Always call copy_to_old here + ! (regardless of use_superad_reduction) so the _old array gets + ! allocated on the first step and check_sizes does not trip on a + ! null pointer downstream. The have_superad_reduction_factor flag + ! gates whether the limiter actually consumes the snapshot. + if (associated(s% superad_reduction_factor)) then + call copy_to_old(s% superad_reduction_factor, & + s% superad_reduction_factor_old, ierr) + if (ierr /= 0) return + if (s% use_superad_reduction) s% have_superad_reduction_factor = .true. + end if + 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) diff --git a/star/private/init.f90 b/star/private/init.f90 index 72b2605fe..421863622 100644 --- a/star/private/init.f90 +++ b/star/private/init.f90 @@ -538,6 +538,7 @@ subroutine set_starting_star_data(s, ierr) s% okay_to_set_mixing_info = .true. s% okay_to_set_mlt_vc = .false. ! not until have set mlt_cv_old s% have_mlt_vc = .false. + s% have_superad_reduction_factor = .false. s% have_ST_start_info = .false. s% prev_mesh_have_ST_start_info = .false. diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90 index 689974af6..3c2e4a0de 100644 --- a/star/private/remove_shells.f90 +++ b/star/private/remove_shells.f90 @@ -1234,6 +1234,7 @@ subroutine do_relax_to_star_cut( & ! save have_mlt_vc and set to false (to load ZAMS model) save_have_mlt_vc = s% have_mlt_vc s% have_mlt_vc = .false. + s% have_superad_reduction_factor = .false. !save composition and entropy profiles xa(:,:) = s% xa(:,k_remove:s% nz) diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 12ae9125a..502c798df 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -221,7 +221,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg ! these are used by use_superad_reduction real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit - real(dp) :: vc_old_local, vc_old_floor + real(dp) :: vc_old_local, vc_old_floor, Gamma_factor_old_local type(auto_diff_real_star_order1) :: tau_conv, f_turnover type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, & diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq, mlt_Pturb, Ptot @@ -522,29 +522,54 @@ subroutine set_superad_reduction() exp_limit = exp(-lambda_limit*(Gamma_factor-1d0)) Gamma_factor = 2d0*(reduction_limit-1d0)*(1d0/(1d0+exp_limit)-0.5d0)+1d0 end if - ! Convective-turnover-time limiter: smoothly suppress the throttle - ! when dt < tau_conv = scale_height / mlt_vc_old. Uses the - ! previous-step converged vc (real(dp), no autoDiff partials) - ! and the step's dt, both frozen as parameters w.r.t. the - ! inner Newton solve. f_turnover -> 1 as dt/tau_conv -> infty - ! (current behavior); f_turnover -> 0 as dt/tau_conv -> 0 - ! (throttle off). Opt-in via superad_reduction_use_turnover_limit. - ! Skip until previous-step mlt_vc has been populated, otherwise - ! s% mlt_vc_old is unassociated and dereferencing it segfaults. + ! Convective-turnover-time limiter: relax Gamma_factor toward + ! its instantaneous (unlimited) value over the local + ! turnover time tau_conv = scale_height / mlt_vc_old. The + ! "starting point" of the relaxation is the previous step's + ! converged Gamma_factor (s% superad_reduction_factor_old(k)), + ! not 1. This is the physically correct first-order ODE + ! + ! d Gamma_fac / dt = (Gamma_fac_inst - Gamma_fac) / tau_conv + ! + ! discretized as + ! + ! Gamma_fac_new = Gamma_fac_old + ! + f_turnover * (Gamma_fac_inst - Gamma_fac_old) + ! + ! When dt >> tau_conv, f_turnover -> 1 and we get the + ! instantaneous value (current behavior). When dt << tau_conv, + ! f_turnover -> 0 and Gamma_fac stays at its previous-step + ! value -- convection had no time to adapt, so neither + ! tightening nor releasing the throttle is allowed. The + ! previous formula relaxed toward 1 (no throttle) in the + ! second limit, which is unphysical -- convection cannot + ! abruptly "release" the throttle either, and that produced + ! HRD jumps when dt suddenly collapsed in late phases. + ! + ! Skip the relaxation until the previous-step mlt_vc and + ! Gamma_factor_old have both been populated. if (s% superad_reduction_use_turnover_limit .and. & Gamma_factor > 1d0 .and. k > 0 .and. & s% have_mlt_vc .and. associated(s% mlt_vc_old) .and. & + s% have_superad_reduction_factor .and. & + associated(s% superad_reduction_factor_old) .and. & s% dt > 0d0) then ! Optional floor on mlt_vc_old at a fraction of the local - ! face sound speed. Prevents tau_conv from blowing up in - ! slow-convection iron-bump cells where f_turnover would - ! otherwise drop to ~0 and leave the throttle off. + ! face sound speed. See controls.defaults docstring for + ! the physical motivation (MLT is unreliable in + ! radiation-pressure-dominated layers; perturbations + ! still propagate at >= c_s). vc_old_floor = s% superad_reduction_turnover_vc_floor_frac & * s% csound_face(k) vc_old_local = max(s% mlt_vc_old(k), vc_old_floor, 1d-30) tau_conv = scale_height / vc_old_local f_turnover = 1d0 - exp(-s% dt / tau_conv) - Gamma_factor = 1d0 + f_turnover * (Gamma_factor - 1d0) + ! Anchor the relaxation at the previous step's value. + ! Clamp to >= 1 (= no-throttle floor) so we never start + ! relaxing from a non-physical sub-1 anchor. + Gamma_factor_old_local = max(s% superad_reduction_factor_old(k), 1d0) + Gamma_factor = Gamma_factor_old_local + & + f_turnover * (Gamma_factor - Gamma_factor_old_local) end if end if end if diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index c69242038..b6d7194a7 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -570,7 +570,14 @@ real(dp), pointer, dimension(:,:) :: xh_old, xa_old real(dp), pointer, dimension(:) :: & - dq_old, q_old, j_rot_old, omega_old, mlt_vc_old + dq_old, q_old, j_rot_old, omega_old, mlt_vc_old, & + superad_reduction_factor_old + ! superad_reduction_factor_old: previous step's converged + ! Gamma_factor (= 1 when no throttle was active). Used by the + ! turnover-time limiter to relax Gamma_factor smoothly toward + ! the new instantaneous target rather than snapping it back to 1 + ! when dt collapses below tau_conv (which would cause an + ! unphysical HRD jump). Companion to mlt_vc_old. real(dp) :: time_old, dt_old, mstar_old, xmstar_old, mstar_dot_old, & M_center_old, v_center_old, R_center_old, L_center_old integer :: nz_old, model_number_old @@ -906,7 +913,11 @@ RSP_just_set_velocities, & have_new_cz_bdy_info, using_gold_tolerances, need_to_setvars, have_new_generation, & using_velocity_time_centering, okay_to_set_mixing_info, & - okay_to_set_mlt_vc, have_mlt_vc + okay_to_set_mlt_vc, have_mlt_vc, & + have_superad_reduction_factor ! true once a step has set + ! superad_reduction_factor and we've + ! copied it into _old (consumed by + ! the smooth-relaxation limiter) real(dp), pointer :: xa_removed(:) ! mass fractions for removed mass (1:species) real(dp) :: h1_czb_mass From dcbb6f3a9d6a062bd98cee0f7bc35ae32fa5a4a7 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Thu, 21 May 2026 19:44:01 -0400 Subject: [PATCH 10/14] docs/changelog: note smooth-relaxation formula for the turnover limiter --- docs/source/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 3c079280e..b0dd49652 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -36,6 +36,8 @@ Opt-in convective-turnover-time limiter for ``use_superad_reduction``: when the 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 From 3c9a152f4332059029fd4aceefd33176526f7842 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Thu, 21 May 2026 19:59:33 -0400 Subject: [PATCH 11/14] docs/controls.defaults: fix sphinx-lint dangling-hyphen at line 2662 CI's sphinx-lint flagged 'density-' at end of line as a dangling hyphen. Move 'density-inversion' onto the next line as a single compound word. Pure docstring fix, no semantic change. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/defaults/controls.defaults | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index c5848219c..971e71d8a 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2659,8 +2659,8 @@ ! 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 + ! 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++. From b3d39da157482bc0475ab44735ef4e49bc557081 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Thu, 21 May 2026 23:43:50 -0400 Subject: [PATCH 12/14] star/superad_reduction: thread Gamma_factor through mesh remap; smooth-relax also on throttle-release; protect snapshot Three correctness bugs in the smooth-relaxation limiter, all flagged in review: 1) Mesh remap. The previous-step Gamma_factor was kept in s%superad_reduction_factor_old and read by the limiter as Gamma_factor_old(k). But s%superad_reduction_factor was not threaded through do_mesh_adjust / do_prune_mesh_surface / prev_mesh restoration alongside mlt_vc. After a split/merge/revised mesh, k on the new mesh no longer mapped to the same mass coordinate as Gamma_factor_old(k) -- the limiter was reading a stale-mesh value. Fixed by mirroring the full mlt_vc plumbing for superad_reduction_factor: - Added prev_mesh_superad_reduction_factor(:) to star_data/public/star_data_step_work.inc. - Allocation in alloc.f90 alongside prev_mesh_mlt_vc. - Save (line 1866-style) and retry-restoration (line 2200-style) in evolve.f90. - Added to do_mesh_adjust signature; interpolated onto the new mesh via do_interp_pt_val (mesh_adjust.f90:~250) with a neutral fill value of 1d0 (= no throttle) for cells without overlap. - Added to do_prune_mesh_surface signature; pruned via prune1. - prev_mesh struct restoration in adjust_mesh.f90:~431. - Threaded through the callers in adjust_mesh.f90 and remove_shells.f90. 2) Throttle-release. The limiter guard was if (... .and. Gamma_factor > 1d0 .and. ...) then meaning the relaxation block only fired when the new step needed throttling instantaneously. If the previous step was throttled (Gamma_factor_old > 1) and the new step had no instantaneous trigger (Gamma_factor = 1), the limiter was skipped and the throttle snapped to 1 in a single step -- exactly the abrupt release the smooth-relaxation form was supposed to prevent. The fix: - Moved the limiter block out of the enclosing `if (Gamma_term > 0d0)` block so it executes regardless of whether the instantaneous trigger fired. - Replaced the Gamma_factor > 1d0 guard with `Gamma_factor > 1d0 .or. Gamma_factor_old > 1d0`, so the relaxation fires whenever there is throttling to evolve -- either tightening it OR releasing it. 3) Start-of-step overwrite. set_vars_if_needed is called at line 1891 of evolve.f90, before new_generation snapshots the previous step's superad_reduction_factor into superad_reduction_factor_old. That path eventually calls Get_results -> set_superad_reduction, which writes the new step's value into s%superad_reduction_factor. So without protection, the snapshot taken at line 1895 would already be the new step's value, not the previous step's converged value. Fix: added s%okay_to_set_superad_reduction_factor flag (mirror of s%okay_to_set_mlt_vc) that gates the write to s%superad_reduction_factor(k) at the end of set_superad_reduction. The flag is set false at all the same sites that set okay_to_set_mlt_vc false (init.f90, evolve.f90:321, evolve.f90:700, remove_shells.f90 when removing shells); set true at the same sites that flip okay_to_set_mlt_vc true (evolve_support.f90:209 right after copy_to_old, read_model.f90:158 after model load). Validation: 60 M_sun to_cc phase restarted from after_core_c_burn completes cleanly to fe_core_infall_limit, no SIGSEGVs on the new arrays, no star_info_old_arrays size-check failures across mesh adjustments. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/private/adjust_mesh.f90 | 2 + star/private/alloc.f90 | 2 + star/private/evolve.f90 | 4 + star/private/evolve_support.f90 | 1 + star/private/init.f90 | 1 + star/private/mesh_adjust.f90 | 19 ++++- star/private/read_model.f90 | 1 + star/private/remove_shells.f90 | 2 + star/private/turb_support.f90 | 99 ++++++++++++------------ star_data/public/star_data_step_work.inc | 12 ++- 10 files changed, 89 insertions(+), 54 deletions(-) diff --git a/star/private/adjust_mesh.f90 b/star/private/adjust_mesh.f90 index a243f2442..1d46dbefd 100644 --- a/star/private/adjust_mesh.f90 +++ b/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) diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 33af7dbef..9c3c58e8c 100644 --- a/star/private/alloc.f90 +++ b/star/private/alloc.f90 @@ -1403,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 diff --git a/star/private/evolve.f90 b/star/private/evolve.f90 index dece07d13..a92b92455 100644 --- a/star/private/evolve.f90 +++ b/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. ! mirror for superad_reduction_factor 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,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 diff --git a/star/private/evolve_support.f90 b/star/private/evolve_support.f90 index f463a22d2..c25be0e82 100644 --- a/star/private/evolve_support.f90 +++ b/star/private/evolve_support.f90 @@ -207,6 +207,7 @@ subroutine set_current_to_old(s) s% mlt_vc(k) = s% mlt_vc_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 diff --git a/star/private/init.f90 b/star/private/init.f90 index 421863622..dbcd3b992 100644 --- a/star/private/init.f90 +++ b/star/private/init.f90 @@ -539,6 +539,7 @@ subroutine set_starting_star_data(s, ierr) s% okay_to_set_mlt_vc = .false. ! not until have set mlt_cv_old s% have_mlt_vc = .false. s% have_superad_reduction_factor = .false. + s% okay_to_set_superad_reduction_factor = .false. s% have_ST_start_info = .false. s% prev_mesh_have_ST_start_info = .false. diff --git a/star/private/mesh_adjust.f90 b/star/private/mesh_adjust.f90 index 3e939771a..c9e0b9ea6 100644 --- a/star/private/mesh_adjust.f90 +++ b/star/private/mesh_adjust.f90 @@ -48,6 +48,7 @@ subroutine do_mesh_adjust( & j_rot_old, i_rot_old, omega_old, D_omega_old, & mlt_vc_old, lnT_old, w_old, specific_PE_old, specific_KE_old, & old_m, old_r, old_rho, dPdr_dRhodr_info_old, D_mix_old, & + superad_reduction_factor_old, & cell_type, comes_from, dq_old, xq_old, xh, xa, dq, xq, ierr) use chem_lib, only: basic_composition_info use interp_1d_def @@ -62,7 +63,8 @@ subroutine do_mesh_adjust( & lnd_old, lnPgas_old, mlt_vc_old, lnT_old, w_old, & specific_PE_old, specific_KE_old, & old_m, old_r, old_rho, dPdr_dRhodr_info_old, & - j_rot_old, omega_old, D_omega_old, D_mix_old + j_rot_old, omega_old, D_omega_old, D_mix_old, & + superad_reduction_factor_old type(auto_diff_real_star_order1), dimension(:), pointer :: i_rot_old real(dp), dimension(:,:), pointer :: xh_old, xa_old real(dp), dimension(:,:), pointer :: xh, xa @@ -247,6 +249,16 @@ subroutine do_mesh_adjust( & 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr) if (failed('mlt_cv')) return + ! Remap the previous-step superadiabatic-reduction factor onto the + ! new mesh so the smooth-relaxation limiter reads the correct local + ! Gamma_factor_old(k) after mesh adjustments. Neutral fill = 1d0 + ! (no throttle) for cells that have no overlap with the old mesh. + call do_interp_pt_val( & + s, nz, nz_old, nzlo, nzhi, & + s% superad_reduction_factor, superad_reduction_factor_old, & + 1d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr) + if (failed('superad_reduction_factor')) return + call do_interp_pt_val( & s, nz, nz_old, nzlo, nzhi, s% D_mix, D_mix_old, & 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr) @@ -556,6 +568,7 @@ subroutine do_prune_mesh_surface( & mlt_vc_old, lnT_old, & dPdr_dRhodr_info_old, nu_ST_old, D_ST_old, D_DSI_old, D_SH_old, & D_SSI_old, D_ES_old, D_GSF_old, D_mix_old, & + superad_reduction_factor_old, & xh, xa, ierr) use auto_diff_support type (star_info), pointer :: s @@ -564,7 +577,8 @@ subroutine do_prune_mesh_surface( & j_rot_old, omega_old, & D_omega_old, am_nu_rot_old, mlt_vc_old, lnT_old, & dPdr_dRhodr_info_old, nu_ST_old, D_ST_old, D_DSI_old, D_SH_old, & - D_SSI_old, D_ES_old, D_GSF_old, D_mix_old + D_SSI_old, D_ES_old, D_GSF_old, D_mix_old, & + superad_reduction_factor_old type(auto_diff_real_star_order1), dimension(:), pointer :: i_rot_old real(dp), dimension(:,:), pointer :: xh_old, xa_old real(dp), dimension(:,:), pointer :: xh, xa @@ -595,6 +609,7 @@ subroutine do_prune_mesh_surface( & call prune1(s% lnT, lnT_old, skip) call prune1(s% D_mix, D_mix_old, skip) call prune1(s% mlt_vc, mlt_vc_old, skip) + call prune1(s% superad_reduction_factor, superad_reduction_factor_old, skip) if (s% rotation_flag) then call prune1(s% j_rot, j_rot_old, skip) diff --git a/star/private/read_model.f90 b/star/private/read_model.f90 index c1967a954..1848457a4 100644 --- a/star/private/read_model.f90 +++ b/star/private/read_model.f90 @@ -157,6 +157,7 @@ subroutine finish_load_model(s, restart, ierr) if (.not. s% have_mlt_vc) then s% okay_to_set_mlt_vc = .true. end if + s% okay_to_set_superad_reduction_factor = .true. s% doing_finish_load_model = .true. call set_vars(s, s% dt, ierr) diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90 index 3c2e4a0de..e702c4b27 100644 --- a/star/private/remove_shells.f90 +++ b/star/private/remove_shells.f90 @@ -1111,6 +1111,7 @@ subroutine do_remove_surface(id, surface_k, ierr) prv% mlt_vc, prv% lnT, & prv% dPdr_dRhodr_info, prv% nu_ST, prv% D_ST, prv% D_DSI, prv% D_SH, & prv% D_SSI, prv% D_ES, prv% D_GSF, prv% D_mix, & + prv% superad_reduction_factor, & s% xh, s% xa, ierr) if (ierr /= 0) then return @@ -1235,6 +1236,7 @@ subroutine do_relax_to_star_cut( & save_have_mlt_vc = s% have_mlt_vc s% have_mlt_vc = .false. s% have_superad_reduction_factor = .false. + s% okay_to_set_superad_reduction_factor = .false. !save composition and entropy profiles xa(:,:) = s% xa(:,k_remove:s% nz) diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 502c798df..f63389663 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -522,59 +522,58 @@ subroutine set_superad_reduction() exp_limit = exp(-lambda_limit*(Gamma_factor-1d0)) Gamma_factor = 2d0*(reduction_limit-1d0)*(1d0/(1d0+exp_limit)-0.5d0)+1d0 end if - ! Convective-turnover-time limiter: relax Gamma_factor toward - ! its instantaneous (unlimited) value over the local - ! turnover time tau_conv = scale_height / mlt_vc_old. The - ! "starting point" of the relaxation is the previous step's - ! converged Gamma_factor (s% superad_reduction_factor_old(k)), - ! not 1. This is the physically correct first-order ODE - ! - ! d Gamma_fac / dt = (Gamma_fac_inst - Gamma_fac) / tau_conv - ! - ! discretized as - ! - ! Gamma_fac_new = Gamma_fac_old - ! + f_turnover * (Gamma_fac_inst - Gamma_fac_old) - ! - ! When dt >> tau_conv, f_turnover -> 1 and we get the - ! instantaneous value (current behavior). When dt << tau_conv, - ! f_turnover -> 0 and Gamma_fac stays at its previous-step - ! value -- convection had no time to adapt, so neither - ! tightening nor releasing the throttle is allowed. The - ! previous formula relaxed toward 1 (no throttle) in the - ! second limit, which is unphysical -- convection cannot - ! abruptly "release" the throttle either, and that produced - ! HRD jumps when dt suddenly collapsed in late phases. - ! - ! Skip the relaxation until the previous-step mlt_vc and - ! Gamma_factor_old have both been populated. - if (s% superad_reduction_use_turnover_limit .and. & - Gamma_factor > 1d0 .and. k > 0 .and. & - s% have_mlt_vc .and. associated(s% mlt_vc_old) .and. & - s% have_superad_reduction_factor .and. & - associated(s% superad_reduction_factor_old) .and. & - s% dt > 0d0) then - ! Optional floor on mlt_vc_old at a fraction of the local - ! face sound speed. See controls.defaults docstring for - ! the physical motivation (MLT is unreliable in - ! radiation-pressure-dominated layers; perturbations - ! still propagate at >= c_s). - vc_old_floor = s% superad_reduction_turnover_vc_floor_frac & - * s% csound_face(k) - vc_old_local = max(s% mlt_vc_old(k), vc_old_floor, 1d-30) - tau_conv = scale_height / vc_old_local - f_turnover = 1d0 - exp(-s% dt / tau_conv) - ! Anchor the relaxation at the previous step's value. - ! Clamp to >= 1 (= no-throttle floor) so we never start - ! relaxing from a non-physical sub-1 anchor. - Gamma_factor_old_local = max(s% superad_reduction_factor_old(k), 1d0) - Gamma_factor = Gamma_factor_old_local + & - f_turnover * (Gamma_factor - Gamma_factor_old_local) - end if end if end if end if - if (k /= 0) s% superad_reduction_factor(k) = Gamma_factor% val + + ! Convective-turnover-time limiter (placed OUTSIDE the + ! Gamma_term>0d0 block so it also handles the case where the + ! instantaneous Gamma_factor is 1 but the previous step's value + ! was > 1 -- i.e., we must relax the throttle DOWN, not snap it + ! to 1 in a single step). Relaxes Gamma_factor toward its + ! instantaneous (capped) value over the local turnover time + ! tau_conv = scale_height / mlt_vc_old, anchored at the + ! previous-step's converged Gamma_factor: + ! + ! d Gamma_fac / dt = (Gamma_fac_inst - Gamma_fac) / tau_conv + ! Gamma_fac_new = Gamma_fac_old + f_turnover * (Gamma_fac_inst - Gamma_fac_old) + ! + ! dt >> tau_conv: f_turnover -> 1, Gamma_fac -> Gamma_fac_inst (no smoothing). + ! dt << tau_conv: f_turnover -> 0, Gamma_fac -> Gamma_fac_old (convection + ! had no time to adapt, throttle held at previous-step value + ! independently of whether the new step needs more or less throttling). + ! + ! Skip until the previous-step mlt_vc and Gamma_factor_old have both + ! been populated. Fires only when there is throttling to track -- + ! either an instantaneous Gamma_factor>1 OR a non-trivial old value + ! that has to relax back to 1. + if (s% superad_reduction_use_turnover_limit .and. k > 0 .and. & + s% have_mlt_vc .and. associated(s% mlt_vc_old) .and. & + s% have_superad_reduction_factor .and. & + associated(s% superad_reduction_factor_old) .and. & + s% dt > 0d0) then + ! Clamp at >= 1 so we never relax from a non-physical sub-1 anchor. + Gamma_factor_old_local = max(s% superad_reduction_factor_old(k), 1d0) + if (Gamma_factor > 1d0 .or. Gamma_factor_old_local > 1d0) then + ! Optional floor on mlt_vc_old at a fraction of the local face + ! sound speed; see controls.defaults for the physical motivation. + vc_old_floor = s% superad_reduction_turnover_vc_floor_frac & + * s% csound_face(k) + vc_old_local = max(s% mlt_vc_old(k), vc_old_floor, 1d-30) + tau_conv = scale_height / vc_old_local + f_turnover = 1d0 - exp(-s% dt / tau_conv) + Gamma_factor = Gamma_factor_old_local + & + f_turnover * (Gamma_factor - Gamma_factor_old_local) + end if + end if + + ! Only commit to s% superad_reduction_factor when allowed; otherwise + ! the start-of-step set_vars_if_needed path would overwrite the + ! previous step's converged value before new_generation has snapshotted + ! it into s% superad_reduction_factor_old. Mirror of okay_to_set_mlt_vc. + if (k /= 0 .and. s% okay_to_set_superad_reduction_factor) then + s% superad_reduction_factor(k) = Gamma_factor% val + end if if (Gamma_factor > 1d0) then grad_scale = (gradr-gradL)/(Gamma_factor*gradr) + gradL/gradr gradr_scaled = grad_scale*gradr diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index b6d7194a7..9540336c9 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -333,7 +333,8 @@ real(dp), pointer, dimension(:,:) :: prev_mesh_xh ! (nvar_hydro,prev_mesh_nz) real(dp), pointer, dimension(:,:) :: prev_mesh_xa ! (species,prev_mesh_nz) real(dp), pointer, dimension(:) :: & ! (prev_mesh_nz) - prev_mesh_j_rot, prev_mesh_omega, prev_mesh_dq, prev_mesh_mlt_vc + prev_mesh_j_rot, prev_mesh_omega, prev_mesh_dq, prev_mesh_mlt_vc, & + prev_mesh_superad_reduction_factor logical :: prev_mesh_species_or_nvar_hydro_changed ! specifies if either species or number of hydro variables ! have been altered since prev_mesh info was stored real(dp), pointer, dimension(:) :: prev_mesh_D_ST_start, prev_mesh_nu_ST_start @@ -914,10 +915,17 @@ have_new_cz_bdy_info, using_gold_tolerances, need_to_setvars, have_new_generation, & using_velocity_time_centering, okay_to_set_mixing_info, & okay_to_set_mlt_vc, have_mlt_vc, & - have_superad_reduction_factor ! true once a step has set + have_superad_reduction_factor, & + ! true once a step has set ! superad_reduction_factor and we've ! copied it into _old (consumed by ! the smooth-relaxation limiter) + okay_to_set_superad_reduction_factor + ! mirror of okay_to_set_mlt_vc: + ! protects s%superad_reduction_factor + ! from being written between the + ! previous step's success and + ! new_generation's snapshot real(dp), pointer :: xa_removed(:) ! mass fractions for removed mass (1:species) real(dp) :: h1_czb_mass From de16d72f90cb22416dab27708cbb380c0c0151f3 Mon Sep 17 00:00:00 2001 From: Matteo Cantiello Date: Fri, 22 May 2026 00:43:39 -0700 Subject: [PATCH 13/14] star/superad_reduction: snapshot Gamma_factor after mesh adjust, drop dead okay_to_set gate The previous commit attempted to protect s%superad_reduction_factor from being overwritten by the start-of-step set_vars_if_needed by mirroring mlt_vc's okay_to_set_mlt_vc gate. That approach failed in practice: the flag is set to .true. only in read_model and in the retry path (set_current_to_old), never in the normal step flow. After step 1, the gate is permanently false and the gated write "s%superad_reduction_factor(k) = Gamma_factor%val" never fires. For mlt_vc this is harmless because s%mlt_vc_ad carries the autoDiff truth and s%mlt_vc is just a plain-dp snapshot; there is no autoDiff back-channel for Gamma_factor, so gating its only write killed it. s%superad_reduction_factor_old ended up frozen at its initial value, the smooth-relaxation anchor (max(_old, 1d0)) collapsed to 1d0, and the formula degenerated back to the old 1+f_tau*(inst-1) behaviour. The HRD endpoint jump returned. Refactored approach: skip the gate entirely and instead snapshot s%superad_reduction_factor -> s%superad_reduction_factor_old at a deterministic point in prepare_for_new_step -- AFTER do_mesh has remapped s%superad_reduction_factor onto the new mesh (so _old is on the correct mesh), and BEFORE set_vars_if_needed could overwrite it with the new step's first-pass value. set s%have_superad_reduction_factor true at the same point. Effect: - The write in set_superad_reduction is now unconditional (matches the original behaviour pre-PR), so s%superad_reduction_factor stays live across Newton iterations. - copy_to_old in new_generation no longer manages the snapshot; it remains in place to keep s%superad_reduction_factor_old allocated (so the existing check_sizes path is happy). - All okay_to_set_superad_reduction_factor sites removed: declaration, init.f90, evolve.f90:321, evolve.f90:702, evolve_support.f90:210, read_model.f90:160, remove_shells.f90:1239. Validation: 60 M_sun full pipeline (60M_cc_smooth_fixed) reaches fe_core_infall_limit cleanly with the relaxation anchored at the real previous-step value. Track on the HRD now differs measurably from the old-formula 60M_cc track, as expected -- the smooth-relaxation is actually doing work this time. Co-Authored-By: Claude Opus 4.7 (1M context) --- star/private/evolve.f90 | 17 +++++++++++++++-- star/private/evolve_support.f90 | 15 ++++++--------- star/private/init.f90 | 1 - star/private/read_model.f90 | 1 - star/private/remove_shells.f90 | 1 - star/private/turb_support.f90 | 8 +------- star_data/public/star_data_step_work.inc | 8 +------- 7 files changed, 23 insertions(+), 28 deletions(-) diff --git a/star/private/evolve.f90 b/star/private/evolve.f90 index a92b92455..f936707ad 100644 --- a/star/private/evolve.f90 +++ b/star/private/evolve.f90 @@ -319,7 +319,6 @@ 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. ! mirror for superad_reduction_factor if (s% timestep_hold > s% model_number + 10000) then write(*,3) 'ERROR: s% timestep_hold', s% timestep_hold, s% model_number @@ -699,7 +698,6 @@ 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 @@ -1891,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 diff --git a/star/private/evolve_support.f90 b/star/private/evolve_support.f90 index c25be0e82..68ffec86f 100644 --- a/star/private/evolve_support.f90 +++ b/star/private/evolve_support.f90 @@ -58,18 +58,16 @@ subroutine new_generation(s, ierr) call copy_to_old(s% mlt_vc, s% mlt_vc_old, ierr) if (ierr /= 0) return - ! Smooth-relaxation limiter for use_superad_reduction needs - ! the previous step's converged Gamma_factor. Mirror the - ! mlt_vc -> mlt_vc_old pattern. Always call copy_to_old here - ! (regardless of use_superad_reduction) so the _old array gets - ! allocated on the first step and check_sizes does not trip on a - ! null pointer downstream. The have_superad_reduction_factor flag - ! gates whether the limiter actually consumes the snapshot. + ! Always ensure superad_reduction_factor_old is allocated, even + ! when use_superad_reduction is .false., so check_sizes does not + ! trip on a null pointer. The actual snapshot used by the + ! smooth-relaxation limiter is done in evolve.f90 between the + ! mesh adjust and set_vars_if_needed (and the + ! have_superad_reduction_factor flag is set there). if (associated(s% superad_reduction_factor)) then call copy_to_old(s% superad_reduction_factor, & s% superad_reduction_factor_old, ierr) if (ierr /= 0) return - if (s% use_superad_reduction) s% have_superad_reduction_factor = .true. end if call enlarge_if_needed_2(s% xh_old,s% nvar_hydro,nz,nz_alloc_extra,ierr) @@ -207,7 +205,6 @@ subroutine set_current_to_old(s) s% mlt_vc(k) = s% mlt_vc_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 diff --git a/star/private/init.f90 b/star/private/init.f90 index dbcd3b992..421863622 100644 --- a/star/private/init.f90 +++ b/star/private/init.f90 @@ -539,7 +539,6 @@ subroutine set_starting_star_data(s, ierr) s% okay_to_set_mlt_vc = .false. ! not until have set mlt_cv_old s% have_mlt_vc = .false. s% have_superad_reduction_factor = .false. - s% okay_to_set_superad_reduction_factor = .false. s% have_ST_start_info = .false. s% prev_mesh_have_ST_start_info = .false. diff --git a/star/private/read_model.f90 b/star/private/read_model.f90 index 1848457a4..c1967a954 100644 --- a/star/private/read_model.f90 +++ b/star/private/read_model.f90 @@ -157,7 +157,6 @@ subroutine finish_load_model(s, restart, ierr) if (.not. s% have_mlt_vc) then s% okay_to_set_mlt_vc = .true. end if - s% okay_to_set_superad_reduction_factor = .true. s% doing_finish_load_model = .true. call set_vars(s, s% dt, ierr) diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90 index e702c4b27..4f404ea28 100644 --- a/star/private/remove_shells.f90 +++ b/star/private/remove_shells.f90 @@ -1236,7 +1236,6 @@ subroutine do_relax_to_star_cut( & save_have_mlt_vc = s% have_mlt_vc s% have_mlt_vc = .false. s% have_superad_reduction_factor = .false. - s% okay_to_set_superad_reduction_factor = .false. !save composition and entropy profiles xa(:,:) = s% xa(:,k_remove:s% nz) diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index f63389663..393a982e5 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -567,13 +567,7 @@ subroutine set_superad_reduction() end if end if - ! Only commit to s% superad_reduction_factor when allowed; otherwise - ! the start-of-step set_vars_if_needed path would overwrite the - ! previous step's converged value before new_generation has snapshotted - ! it into s% superad_reduction_factor_old. Mirror of okay_to_set_mlt_vc. - if (k /= 0 .and. s% okay_to_set_superad_reduction_factor) then - s% superad_reduction_factor(k) = Gamma_factor% val - end if + if (k /= 0) s% superad_reduction_factor(k) = Gamma_factor% val if (Gamma_factor > 1d0) then grad_scale = (gradr-gradL)/(Gamma_factor*gradr) + gradL/gradr gradr_scaled = grad_scale*gradr diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 9540336c9..f355021a8 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -915,17 +915,11 @@ have_new_cz_bdy_info, using_gold_tolerances, need_to_setvars, have_new_generation, & using_velocity_time_centering, okay_to_set_mixing_info, & okay_to_set_mlt_vc, have_mlt_vc, & - have_superad_reduction_factor, & + have_superad_reduction_factor ! true once a step has set ! superad_reduction_factor and we've ! copied it into _old (consumed by ! the smooth-relaxation limiter) - okay_to_set_superad_reduction_factor - ! mirror of okay_to_set_mlt_vc: - ! protects s%superad_reduction_factor - ! from being written between the - ! previous step's success and - ! new_generation's snapshot real(dp), pointer :: xa_removed(:) ! mass fractions for removed mass (1:species) real(dp) :: h1_czb_mass From d88fe3c41800357447a5ce8e9c9b68760855c0df Mon Sep 17 00:00:00 2001 From: Matthias Fabry Date: Tue, 2 Jun 2026 00:40:59 -0400 Subject: [PATCH 14/14] rebase to Matteo's branch, add optional linear vs exponential (#1008) --- star/defaults/controls.defaults | 20 +++++++++++++++----- star/private/ctrls_io.f90 | 5 ++++- star/private/turb_support.f90 | 7 ++++++- star_data/private/star_controls.inc | 1 + 4 files changed, 26 insertions(+), 7 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 971e71d8a..950aec433 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2746,11 +2746,8 @@ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Optional convective-turnover-time limiter for the throttle. - ! When ``.true.``, multiplies the throttle excess - ! ``(Gamma_factor - 1)`` by - ! - ! ``f_tau = 1 - exp(-dt / tau_conv)`` - ! + ! 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 @@ -2776,6 +2773,19 @@ 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 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index e12508ed1..507a4bd0e 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -136,7 +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_vc_floor_frac, & + 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, & @@ -1080,6 +1081,7 @@ subroutine store_controls(s, ierr) 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 @@ -2805,6 +2807,7 @@ subroutine set_controls_for_writing(s, ierr) 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 diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 393a982e5..834d5ceec 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -561,7 +561,12 @@ subroutine set_superad_reduction() * s% csound_face(k) vc_old_local = max(s% mlt_vc_old(k), vc_old_floor, 1d-30) tau_conv = scale_height / vc_old_local - f_turnover = 1d0 - exp(-s% dt / tau_conv) + select case (trim(s% superad_reduction_turnover_limit_function)) + case ('linear') + f_turnover = min(s% dt / tau_conv, 1d0) + case default + f_turnover = 1d0 - exp(-s% dt / tau_conv) + end select Gamma_factor = Gamma_factor_old_local + & f_turnover * (Gamma_factor - Gamma_factor_old_local) end if diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index bd86209b1..56e8b19bd 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -307,6 +307,7 @@ superad_reduction_diff_grads_limit, & superad_reduction_limit logical :: superad_reduction_use_turnover_limit + character(len=strlen) :: superad_reduction_turnover_limit_function real(dp) :: superad_reduction_turnover_vc_floor_frac ! mixing parameters