superad turnover time limiter - #1045
Open
Debraheem wants to merge 20 commits into
Open
Conversation
…ibution 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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
… 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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
…_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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
…ep 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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
…h-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) <noreply@anthropic.com>
… 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) <noreply@anthropic.com>
…rnover limiter The turnover limiter interpolated Gamma_factor toward its instantaneous value. But the applied throttle is eta = 1/Gamma_factor (gradr_scaled = gradL + eta*(gradr-gradL)), so the suppression is linear in eta, not Gamma. Relax eta and invert back, so f_turnover is the fraction of the *applied* reduction per step. eta stays in (0,1] (anchor clamped to Gamma_old >= 1), preserving Gamma_factor >= 1. Composes with the linear/exponential option from #1008. Validation: 60 Msun ZAMS->core-He benchmarks (limiter on, +/- v_c floor) reach the same endpoints with no regression -- retries and throttle within ~0.5% of the old behavior. Those runs keep the throttle saturated near the cap, so the change is felt mainly on low->high throttle transitions. Thanks to @Debraheem for catching this in review.
Debraheem
changed the title
Member
Author
|
I made this pr separate because i did not want to overwrite the previous pr opened by Matteo until this branch is well tested. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This pr attempts to clean up and refactor #1006. Still needs further testing and could be polished and furthered cleaned.
For now I replaced the mlt++ in the 80M test case with this. https://testhub.mesastar.org/EbF_matteo_matthias%2Fsuperad_turnover_time_limiter/commits/a5c535a