getfc.f

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

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

Variables on this page (6)

cc13, cln, dm, fc, json, rm

Purpose

GETFC computes, shell by shell, the factor $f_c$ that relates the diffusion coefficient for chemical mixing to that for angular momentum transport by meridional circulation, "using the prescription of Zahn (1992)" (its comment).

How it fits

Who calls it: only SECULR, when LVFC is set, for the shells IMIN..IMAX.
Inputs: HD ($\log\rho$ at zone centres), HRU (radii at zone centres), HV (circulation velocities $U$ at zone edges), OMEGA; and from INTVAR_MOD, RM (radii at zone edges, the means of neighbouring centre radii) and DM (density at zone edges, interpolated in $\log\rho$), both set by SETUPV.
Output: the argument VFC, an array that SECULR owns (it clears it, and passes it on to CODIFF, which reads it).
Modules: NMLPAR_MOD gives FC and LVFC (namelist parameters, read in PARMIN); INTVAR_MOD (the edge quantities), shared with CODIFF, SECULR, SETUPV, VCIRC; CONSTANTS_MOD (CLN, CC13); DIMENSIONS_MOD (JSON).

The method

The comment gives $$f_c = \min\left(1,\ \frac{C}{30}\,\frac{U}{2V - \alpha U}\right),\qquad \alpha = \frac12\frac{d\ln(r^2\Omega)}{d\ln r},\qquad V = \frac{1}{6\rho r}\frac{d}{dr}(\rho r^2 U).$$ The code computes:

ALPHA $= \dfrac{1}{\Omega_m r_m}\dfrac{\Delta(\Omega r^2)}{\Delta r}$ between neighbouring zone centres, with $\Omega_m$ the mean of the two and $r_m$ = RM the edge between them, which is $d\ln(r^2\Omega)/d\ln r$, without the factor ½ (Bugs, by priority, 1). It uses shell IMIN-1, which exists: MIDMOD sets IMIN ≥ 2.
V $= \dfrac{1}{3\rho r}\dfrac{d}{dr}(\rho r^2 U)$, which is $2V$; at the first and last shells, by the product rule with $dU/dr$ from the neighbour.
VFC $= \dfrac{f_c^{\mathrm{in}}}{\max\left(|{\tt V} - {\tt ALPHA}\,U|/U,\ f_c^{\mathrm{in}}\right)} = \min\left(1,\ \dfrac{f_c^{\mathrm{in}}\,U}{|2V - {\tt ALPHA}\,U|}\right)$, with $f_c^{\mathrm{in}}$ = FC; zero where $U \le 0$, and everywhere if FC ≤ 0. The absolute value is the code's; the comment's formula has none.

Question 1 (line 63 ↗, lines 112–113 ↗). FC and $C/30$. The input FC plays the part of $C/30$ in the comment's formula, while without LVFC it is the constant ratio of the two coefficients (CODIFF). Is that the intended meaning of the parameter with LVFC?

Question 2 (lines 109–110 ↗). Negative velocities. Where $U \le 0$, $f_c = 0$: no chemical mixing where the circulation flows inward. Intended, or should $|U|$ be used?

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. (changes results, small) Checked by KC (Nov 10, 2026) The factor ½ in $\alpha$.

(line 63 ↗; lines 75–79 ↗; line 112 ↗) The comment defines $\alpha = \tfrac12\,d\ln(J/M)/d\ln r$, as in Zahn (1992), and uses $2V - \alpha U$. The code computes ALPHA $= d\ln(r^2\Omega)/d\ln r$, with no ½, and then V - ALPHA*HV, where V is already $2V$. So the code's denominator is $2V - 2\alpha U$, and the $\alpha$ term is twice what the comment gives, which changes $f_c$ wherever that term matters. Only with LVFC (off by default; "experimental" in the Physics namelist page).

A possible fix (unless the comment is the slip):

getfc.f

@@ lines 72–80 @@
          GOTO 9999
       ENDIF
 C DETERMINE ALPHA.
+C 11/09/26 KC via Claude
+C code changed to include the factor 1/2 of the comment's
+C ALPHA=1/2 D(LN J/M)/D LN R.
       DO I = IMIN,IMAX
          WM = 0.5D0*(OMEGA(I)+OMEGA(I-1))
-         ALPHA(I) = (OMEGA(I)*HRU(I)**2-OMEGA(I-1)*HRU(I-1)**2)/
+         ALPHA(I) = 0.5D0*(OMEGA(I)*HRU(I)**2-OMEGA(I-1)*HRU(I-1)**2)/
      *              (HRU(I)-HRU(I-1))/(WM*RM(I))
       END DO
 C DETERMINE 2V.

How the code itself shows it. The routine carries the comment's factors into V but not into ALPHA. The comment defines $V = 1/(6\rho r)\,d(\rho r^2 U)/dr$ (line 64 ↗); the loop under "DETERMINE 2V." (line 80 ↗) computes V(I) = CC13*(DR2VP-DR2VM)/RHOR/(HRU(I)-HRU(I-1)) (line 89 ↗), with CC13 = 1/3, twice the comment's 1/6, and line 99 ↗ and line 105 ↗ do the same at the end shells. For ALPHA, lines 77–78 ↗ compute $d(\Omega r^2)/dr$ divided by $\Omega r$, that is $d\ln(r^2\Omega)/d\ln r$. YREC's $J/M$ is $\Omega$ times $I/M$ (HJM(J) = OMEGA(J)*HI(J), getnewenv.f, line 297 ↗), with $I/M = \tfrac23 r^2$ as MIDMOD's first guess (HIM(I) = CC23*EXP(CLN*2.0D0*HRM(I)), midmod.f, line 493 ↗), so ALPHA is the comment's $d\ln(J/M)/d\ln r$ in full, without the ½ of line 63 ↗. line 112 ↗ then takes V(I)-ALPHA(I)*HV(I): the comment's $2V$ with twice its $\alpha U$. Is the comment's ½, the factor of Zahn (1992) it cites, the one intended?

Checked by running YREC (Nov 9 and 10, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): None of the five rotating test cases (testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2, testsuite/Test_solarGS98_difrotmixper18p4d1p0, testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, testsuite/Test_m2p0_feh+0p0_vzams150_IES1_TAMS, and testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS) reaches it, since none sets LVFC. No namelist in the repository does: Those that list it set it false, with the comment "Not recommended (disables FC setting.)," and the documentation lists it among the experimental options (docs/source/namelist_physics.md, line 298). So on Nov 10, testsuite/Test_m2p0_feh+0p0_vzams150_IES1_TAMS (2 $M_\odot$, solid-body rotation at about 148 km/s at the start, run from the ZAMS to the TAMS) was run with only LVFC switched on: a documented option that YREC's own namelists call "Not recommended," which replaces the case's constant FC = 0.033 with the variable factor computed here. GETFC then computed $f_c$ at 536,252 shells over the run (shells 2 to 758), at every iteration of SECULR (seculr.f, line 192 ↗, seculr.f, line 232 ↗). At all but 370 of them, the factor ½ of the comment would change $f_c$ by more than 1%, by factors from 0.03 to 31. With the factor ½ (the fix above), the star mixes more: The main sequence lasts 1.357 Gyr instead of 1.344 Gyr (635 models instead of 624); before central $X$ falls below 0.05, surface $^{14}$N rises by up to 0.0039 dex, against $7 \times 10^{-6}$ dex for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$, the noise floor); surface $X$ falls by up to $4.9 \times 10^{-5}$, against $4.8 \times 10^{-7}$; central $X$ differs by up to 0.018, against 0.002; and $\log L$ by up to 0.0075 dex, against 0.0011 dex. The effect is small but well above the noise, in a setting YREC does not recommend. Intent: The comment (line 63 ↗) gives the factor ½, so either the code or the comment is the slip.

Definitions

subroutineGETFC

lines 4–119 · GitHub

Zahn (1992) ratio $f_c$ of the chemical to the angular-momentum diffusion coefficient, per shell, capped at 1 (equations above).

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, INTVAR_MOD, NMLPAR_MOD

Called by SECULR ×2

Change log

The changes since rp_1

Changed since rp_1. VFC is a new argument (it was in /VARFC/, now removed; SECULR owns the array); COMMON blocks to modules (FC from NMLPAR_MOD, /VMULT/ gone), IMPLICIT NONE. No question is settled.

1C 2C$$$$$$GETFC 3C 4 SUBROUTINE GETFC(HD,HRU,HV,IMIN,IMAX,OMEGA,VFC) 5C 10/15/26 MHP via Claude 6C code changed to take JSON from DIMENSIONS_MOD instead of declaring 7C PARAMETER(JSON=5000) here. 8 USE DIMENSIONS_MOD 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/INTVAR/ into INTVAR_MOD. 11 USE INTVAR_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 14 USE NMLPAR_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 17 USE CONSTANTS_MOD 18C 10/20/26 MHP via Claude 19C code changed to replace implicit typing with IMPLICIT NONE and 20C explicit declarations. 21 IMPLICIT NONE 22 REAL*8 ALPHA,DR2VM,DR2VP,DRHOR2,DUDR,FTEST,FX,HD,HRU,HV,OMEGA, 23 * * RHOR,V,VFC,VM,VP,WM 24 INTEGER I,IMAX,IMIN 25 26C 10/13/26 MHP via Claude 27C VFC is now a dummy argument, formerly COMMON/VARFC/. 28 DIMENSION VFC(JSON) 29C 9/25/26 MHP via Claude 30C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 31C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 32C 10/12/26 MHP via Claude 33C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 34C namelist variables that the physics routines read from 27 COMMON 35C blocks which also held derived or working members. Those blocks keep 36C only their non-namelist members; five of them -- NEWMX, VMULT, 37C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 38C Derivation-only namelist variables from the same blocks went to 39C COMMON/NMLBYPASS/ instead, not here. 40C 41C Declared identically in all 121 files that hold any of these members. 42C Generated from one member list, so every copy is byte-identical; the 43C linker's 'size of symbol changed' warning is the check that they are. 44C 10/13/26 MHP via Claude 45C code changed to replace COMMON/VARFC/ with an explicit argument: seculr.f 46C clears VFC, GETFC fills it and CODIFF reads it, all inside seculr.f's 47C call tree, so seculr.f now owns the storage. 48 DIMENSION HD(JSON),HRU(JSON),HV(JSON), 49 * * OMEGA(JSON),ALPHA(JSON),V(JSON) 50 SAVE 51C INPUT VARIABLES : 52C HD - LOG10 OF DENSITY(GM/CM**3) AT ZONE CENTERS. 53C RM - UNLOGGED RUN OF RADII AT ZONE EDGES. 54C HRU - UNLOGGED RUN OF RADII AT ZONE CENTERS. 55C HV - UNLOGGED RUN OF DIFFUSION VELOCITIES AT ZONE EDGES. 56C IMIN,IMAX - STARTING AND ENDING SHELLS. 57C OMEGA - RUN OF ANGULAR VELOCITIES AT ZONE CENTERS. 58C OUTPUT VARIABLE : 59C VFC - RATIO OF DIFFUSION COEFFICIENT FOR MIXING TO THAT FOR 60C ANGULAR MOMENTUM TRANSPORT. 61C THIS SUBROUTINE SOLVES FOR THE VARIABLE FC USING THE 62C PRESCRIPTION OF ZAHN (1992). 63C FC = MIN (1, C/30*U/(2V-ALPHA*U), WHERE ALPHA=1/2 D(LN J/M)/D LN R 64C AND V = 1/(6*RHO*R) *D/DR(RHO* R**2 * U) (U = DIFFUSION VELOCITY). 65C 66C IF FC = 0, SKIP CALCULATION (NO MIXING). 67 68 IF(FC.LE.0.0D0)THEN 69 DO I = IMIN,IMAX 70 VFC(I) = 0.0D0 71 END DO 72 GOTO 9999 73 ENDIF 74C DETERMINE ALPHA. 75 DO I = IMIN,IMAX 76 WM = 0.5D0*(OMEGA(I)+OMEGA(I-1)) 77 ALPHA(I) = (OMEGA(I)*HRU(I)**2-OMEGA(I-1)*HRU(I-1)**2)/ 78 * * (HRU(I)-HRU(I-1))/(WM*RM(I)) 79 END DO 80C DETERMINE 2V. 81 FX = (HRU(IMIN)-RM(IMIN))/(RM(IMIN+1)-RM(IMIN)) 82 VM = HV(IMIN)+FX*(HV(IMIN+1)-HV(IMIN)) 83 DR2VM = EXP(CLN*HD(IMIN))*HRU(IMIN)**2*VM 84 DO I = IMIN + 1, IMAX - 1 85 FX = (HRU(I)-RM(I))/(RM(I+1)-RM(I)) 86 VP = HV(I)+FX*(HV(I+1)-HV(I)) 87 DR2VP = EXP(CLN*HD(I))*HRU(I)**2*VP 88 RHOR = DM(I)*RM(I) 89 V(I) = CC13*(DR2VP-DR2VM)/RHOR/(HRU(I)-HRU(I-1)) 90 DR2VM = DR2VP 91 END DO 92C FOR FIRST AND LAST SHELLS, CALCULATE U*D/DR(RHO*R**2) AND ASSUME 93C DU/DR = VALUE DERIVED FROM THE LAST NEIGHBOR SHELL WITH NONZERO V. 94 RHOR = DM(IMIN)*RM(IMIN) 95 DRHOR2 = (EXP(CLN*HD(IMIN))*HRU(IMIN)**2 - 96 * * EXP(CLN*HD(IMIN-1))*HRU(IMIN-1)**2)/ 97 * * (HRU(IMIN) - HRU(IMIN-1)) 98 DUDR = (HV(IMIN+1)-HV(IMIN))/(RM(IMIN+1)-RM(IMIN)) 99 V(IMIN) = CC13*(HV(IMIN)*DRHOR2/RHOR + RM(IMIN)*DUDR) 100 RHOR = DM(IMAX)*RM(IMAX) 101 DRHOR2 = (EXP(CLN*HD(IMAX))*HRU(IMAX)**2 - 102 * * EXP(CLN*HD(IMAX-1))*HRU(IMAX-1)**2)/ 103 * * (HRU(IMAX) - HRU(IMAX-1)) 104 DUDR = (HV(IMAX)-HV(IMAX-1))/(RM(IMAX)-RM(IMAX-1)) 105 V(IMAX) = CC13*(HV(IMAX)*DRHOR2/RHOR + RM(IMAX)*DUDR) 106C NOW COMPUTE RUN OF FC; THIS ASSUMES THAT FC = SAME FOR ALL 107C MECHANISMS AND IS LIMITED TO A MAXIMUM OF 1. 108 DO I = IMIN,IMAX 109 IF(HV(I).LE.0.0D0)THEN 110 VFC(I) = 0.0D0 111 ELSE 112 FTEST = MAX(ABS(V(I)-ALPHA(I)*HV(I))/HV(I),FC) 113 VFC(I) = FC/FTEST 114 ENDIF 115 END DO 116 9999 9999 CONTINUE 117 118 RETURN 119 END

1"""Python version of getfc.f (YREC development, 7be4133): the factor f_c 2relating the chemical diffusion coefficient to that of angular momentum 3transport by meridional circulation (Zahn 1992).""" 4 5import math 6 7import numpy as np 8 9from pyrec import constants_mod 10from pyrec import dimensions_mod 11from pyrec import intvar_mod 12from pyrec import nmlpar_mod 13 14# Locals kept between calls (SAVE). 15alpha = np.zeros(dimensions_mod.json + 1) 16v = np.zeros(dimensions_mod.json + 1) 17 18 19def getfc(hd: np.ndarray, hru: np.ndarray, hv: np.ndarray, imin: int, imax: int, omega: np.ndarray, 20 vfc: np.ndarray) -> None: 21 """Compute, for the shells imin..imax, Zahn's (1992) ratio f_c of the 22 chemical to the angular-momentum diffusion coefficient, capped at 1 23 (equations on the page); zero where the circulation velocity is <= 0, 24 and everywhere if FC <= 0. 25 26 Translates: getfc.f#GETFC (7be4133) 27 28 Arguments: hd (log rho at zone centres), hru (radii at zone centres), 29 hv (circulation velocities U at zone edges), imin, imax, omega. Array 30 changed in place: vfc (owned by seculr, read by codiff). Nothing is 31 returned. 32 33 """ 34 35 c = constants_mod 36 iv = intvar_mod 37 fc = nmlpar_mod.fc 38 39 # No mixing if FC <= 0 (Fortran lines 66–73). 40 if fc <= 0.0: 41 for i in range(imin, imax + 1): 42 vfc[i] = 0.0 43 return # go to 9999 44 45 # alpha = d ln(r^2 Omega) / d ln r between neighbouring centres. 46 # Suspected bug (see the page): no factor 1/2 (Fortran lines 74–79). 47 for i in range(imin, imax + 1): 48 wm = 0.5 * (omega[i] + omega[i - 1]) 49 alpha[i] = ((omega[i] * hru[i]**2 - omega[i - 1] * hru[i - 1]**2) 50 / (hru[i] - hru[i - 1]) / (wm * iv.rm[i])) 51 52 # v = (1/(3 rho r)) d(rho r^2 U)/dr, i.e. 2V; U interpolated to the 53 # zone centres (Fortran lines 80–91). 54 fx = (hru[imin] - iv.rm[imin]) / (iv.rm[imin + 1] - iv.rm[imin]) 55 vm = hv[imin] + fx * (hv[imin + 1] - hv[imin]) 56 dr2vm = math.exp(c.cln * hd[imin]) * hru[imin]**2 * vm 57 for i in range(imin + 1, imax): 58 fx = (hru[i] - iv.rm[i]) / (iv.rm[i + 1] - iv.rm[i]) 59 vp = hv[i] + fx * (hv[i + 1] - hv[i]) 60 dr2vp = math.exp(c.cln * hd[i]) * hru[i]**2 * vp 61 rhor = iv.dm[i] * iv.rm[i] 62 v[i] = c.cc13 * (dr2vp - dr2vm) / rhor / (hru[i] - hru[i - 1]) 63 dr2vm = dr2vp 64 65 # The first and last shells, by the product rule (Fortran lines 66 # 92–105). 67 rhor = iv.dm[imin] * iv.rm[imin] 68 drhor2 = ((math.exp(c.cln * hd[imin]) * hru[imin]**2 69 - math.exp(c.cln * hd[imin - 1]) * hru[imin - 1]**2) 70 / (hru[imin] - hru[imin - 1])) 71 dudr = (hv[imin + 1] - hv[imin]) / (iv.rm[imin + 1] - iv.rm[imin]) 72 v[imin] = c.cc13 * (hv[imin] * drhor2 / rhor + iv.rm[imin] * dudr) 73 rhor = iv.dm[imax] * iv.rm[imax] 74 drhor2 = ((math.exp(c.cln * hd[imax]) * hru[imax]**2 75 - math.exp(c.cln * hd[imax - 1]) * hru[imax - 1]**2) 76 / (hru[imax] - hru[imax - 1])) 77 dudr = (hv[imax] - hv[imax - 1]) / (iv.rm[imax] - iv.rm[imax - 1]) 78 v[imax] = c.cc13 * (hv[imax] * drhor2 / rhor + iv.rm[imax] * dudr) 79 80 # f_c = FC / max(|v - alpha U| / U, FC) (Fortran lines 106–115). 81 for i in range(imin, imax + 1): 82 if hv[i] <= 0.0: 83 # Question (see the page): no mixing where U <= 0. 84 vfc[i] = 0.0 85 else: 86 ftest = max(abs(v[i] - alpha[i] * hv[i]) / hv[i], fc) 87 vfc[i] = fc / ftest