PyREC scaffolding (from YREC development, 7be4133)/ setup_grsett.fGlossary

setup_grsett.f

330 lines · 1 definitions · on GitHub (7be4133)

Python version: setup_grsett.py, 281 lines, on the right (switch to Python above the code).

Variables on this page (19)

cc13, cln, con_mass, con_rad, con_temp, con_time, cstdiffmix, cstmixing, fgry, fgrz, ilambda, ishort, json, ldifz, lthoul, lthoulfit, sdel, xmin, ymin

Purpose

Prepares the gravitational settling (microscopic diffusion) of helium, and optionally of metals, for the older settling solver: it finds the range of shells where settling is computed, converts the model to the units of Bahcall and Loeb, and computes the diffusion coefficients and their derivatives with respect to $X$.

How it fits

  • Who calls it: GRSETT, the Bahcall–Loeb settling solver (used when LNEWDIF is false; see ROTMIX and MIX). The newer path is MICRODIFF_SETUP.
  • Inputs: the step DT (s), $d\ln P/dr$ (HQPR), $\log r$, $\log\rho$, $\log T$ (HR, HD, HT), the enclosed masses HS1 (g), the total mass STOT, the convection flags LC and HCOMP; from modules, SDEL(2,·) (SCRTCH_MOD, which ROTMIX loads with the gradient DELM) and the namelist parameters in NMLPAR_MOD: ILAMBDA, XMIN, YMIN, LTHOULFIT, FGRY, FGRZ, LTHOUL, LDIFZ, CSTMIXING and CSTDIFFMIX.
  • What it calls: THDIFF, the full diffusion coefficients of Thoul et al., when LTHOUL is true and LTHOULFIT false.
  • Outputs: IBEGIN, IEND (the range), LALLCZ (true when settling is skipped), HRU, HTU (radius and temperature in Bahcall–Loeb units), HQPR, HS1, DT, STOT converted in place, and the coefficients COD1, COD2, QCOD1X, QCOD2X; for metals, COD1Z, COD2Z, QCOD1Z, QCOD2Z, arguments whose storage GRSETT owns and passes on to MODEL_TO_EQUAL. The conversion factors go to CONFAC_MOD (CON_RAD, CON_MASS, CON_TEMP, CON_TIME), shared with EQUAL_TO_MODEL and the MICRODIFF routines.

The method

  1. Range. IBEGIN is the first shell above a convective core with $X >$ XMIN; IEND is the last shell below a convective envelope. The helium check then moves IEND down to the first shell at or below it with $Y >$ YMIN (settled from the code, Nov 9, 2026: IEND = I after label 47, lines 149–160↗). A fully convective model, or one with no shell above either limit, is skipped (LALLCZ), with a message; the comment before the helium check's exit says "HYDROGEN-FREE MODEL," but the test and the printed message are on helium (lines 154–156↗). The coefficients below are computed for every shell, 1 to M, not only within the range.

  2. Units (Bahcall and Loeb): $r$ in $R_\odot$ (6.9598×10$^{10}$ cm), $T$ in $10^7$ K, masses multiplied by CON_MASS $=10^{-2}/R_\odot^3$ ($R_\odot$ in cm), time in $2.7\times10^{13}$ yr ("includes factor of 2.2 from ln lambda").

From YREC's documentation

Physics Options (Gravitational settling & microscopic diffusion): LDIFY, LDIFZ, LDIFLI turn on the settling of helium, metals and the light species; FGRY, FGRZ (default 1) scale the diffusion coefficients of helium and metals (for example 0.75 without rotational mixing); XMIN and YMIN (1.0D-3) turn diffusion off at low central $X$ (and, the documentation says, low surface $Y$). ILAMBDA chooses the Coulomb logarithm of the collision integrals: the documentation's line reads "4 (recommended) is a constant, 2.2." Among the rarely used options: LTHOUL, the formalism of Thoul et al., and LTHOULFIT, its published fit instead of the full solution (not recommended).

Suspected bugs, checked by running the code

Each bug is marked by what running the code showed: (changes results, large) or (small), the fix moving a run clearly beyond a null change; (stops the run); (wrong output), the fix changing only what is written; (harmless in the runs), within the noise; (never occurred), the code ran but the fault did not; (not reached); (cannot be reached), a guard prevents it; (withdrawn), the suspicion was wrong.

1. (not reached) The derivatives miss the mixing factors.

(lines 288–291↗; line 311↗; line 316↗) With LTHOUL, COD1 and COD2 are multiplied by CSTMIXING, but their $X$-derivatives QCOD1X, QCOD2X are not; COD1Z is multiplied by CSTDIFFMIX*CSTMIXING, but its derivative QCOD1Z by CSTMIXING only. GRSETT uses the derivatives to correct the coefficients for the change in $X$ (and $Z$) within the step (grsett.f, lines 173–185↗; get_imp_diffco.f, lines 19–21↗), so with a factor other than 1 the correction is that of the unscaled coefficient. What follows: The second-order correction of the settling coefficients is off by the factor; the size of the effect was not estimated. Only with LTHOUL (off by default; forced on by LDIFZ and LDIFLI) and CSTMIXING or CSTDIFFMIX set to a value other than 1 (both default 1).

A possible fix:

setup_grsett.f

@@ lines 285–294 @@
             QACX = -1.42D0 + 1.294D0*X
 CFD 10/09 Mimic Mixing to reduce settling.
 C            COD1(I) = FAC*HQPR(I)*X*(AP+AT)
+C 11/09/26 KC via Claude
+C code changed to scale the X-derivatives by the same mixing factors
+C as the coefficients.
             COD1(I) = CSTMIXING*FAC*HQPR(I)*X*(AP+AT)
             COD2(I) = CSTMIXING*FAC*AC
-            QCOD1X(I) = FAC*HQPR(I)*(AP+AT+X*(QAPX+QATX))
-            QCOD2X(I) = FAC*QACX
+            QCOD1X(I) = CSTMIXING*FAC*HQPR(I)*(AP+AT+X*(QAPX+QATX))
+            QCOD2X(I) = CSTMIXING*FAC*QACX
          ENDIF
 C        METAL DIFFUSION, USING THE THOUL ET AL. COEFFICIENTS FOR FULLY
 C        IONIZED IRON.
@@ lines 313–319 @@
 C old ver.     COD2Z(I) = ABS(FAC*AH)
 C old ver.     QCOD1Z(I) = FAC*HQPR(I)*(AP+AT)
                COD2Z(I) = CSTMIXING*ABS(FAC*AH)
-               QCOD1Z(I) = CSTMIXING*FAC*HQPR(I)*(AP+AT)
+               QCOD1Z(I) = CSTDIFFMIX*CSTMIXING*FAC*HQPR(I)*(AP+AT)
                QCOD2Z(I) = 0.0D0
             ENDIF
          ENDIF

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): not reached in any of 34 runs (26 of YREC's cases, from 0.04 to 9 $M_\odot$, with rotation, settling, and helium burning, and 8 made-up ones). These lines belong to the old settling routines (LNEWDIF false), which no case uses: Every case with settling sets LNEWDIF and LDIFLI, and the documentation lists LNEWDIF among the legacy options (docs/source/namelist_physics.md, line 285). Reaching them takes LDIFY or LDIFZ with LNEWDIF false and LDIFLI off, since LDIFLI forces LNEWDIF on (parmin.f, lines 1239–1251↗), with LTHOUL and CSTMIXING or CSTDIFFMIX other than 1; no case sets either, and the documentation does not mention them. Intent: The comments "CFD 10/09 Mimic Mixing to reduce settling" (line 286↗) and "old ver. QCOD1Z(I) = FACHQPR(I)(AP+AT)" (line 314↗) show that the coefficients and QCOD1Z were edited to take the factors and QCOD1X and QCOD2X were not, which looks like an unfinished edit.

Definitions

subroutineSETUP_GRSETT

lines 2–330 · GitHub

Uses modules CONFAC_MOD, CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD, SCRTCH_MOD

Called by GRSETT

Flowchart of its calls, in order (1)
for I = 1,M (DO 60), lines 189–327
else, lines 223–292
if (.NOT.LTHOULFIT), lines 227–277
THDIFF275

Change log

The changes since rp_1

Changed since rp_1. The metals' coefficients COD1Z, COD2Z, QCOD1Z, QCOD2Z are arguments (formerly /GRAVSZ/); the unused /GSCOF/ is dropped; otherwise declarations only: COMMON blocks to modules, IMPLICIT NONE.