equal_to_model.f

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

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

Variables on this page (7)

con_mass, con_rad, con_temp, con_time, ez, json, ldifz

Purpose

EQUAL_TO_MODEL ends a step of gravitational settling (microscopic diffusion): the diffusion equations are solved on an equally spaced grid in radius, and this routine applies the resulting changes in hydrogen (and, with metal diffusion, in $Z$) to the model's own points. It then converts the model's arrays back from the dimensionless units used by the diffusion solver.

How it fits

Who calls it: only GRSETT, at the end of the settling step; MODEL_TO_EQUAL does the opposite transformation at its start. MICRODIFF_ETM has the same unit conversion for the newer diffusion scheme.
Inputs: the equally spaced radii ER (NPT points) and the changes EX in $X$ on them; the range IBEGIN..IEND of model points inside the grid; the model arrays HCOMP, HRU (radius), HTU (temperature), HS1 (mass), HQPR, the time step DT and the total mass STOT, all in the solver's units.
What it calls: INTRP2, four-point Lagrange weights.
Outputs: HCOMP updated, and HRU, HTU, HS1, HQPR, DT, STOT back in physical units.
Modules: CONFAC_MOD holds the conversion factors CON_RAD, CON_MASS, CON_TEMP, CON_TIME, set by MICRODIFF_SETUP and SETUP_GRSETT; GRAVEZ_MOD holds the change in $Z$ on the grid, EZ (with GRSETT, MODEL_TO_EQUAL and the Lax–Wendroff routines, which use its other entries; this routine reads only EZ); NMLPAR_MOD gives the flag LDIFZ (formerly in /GRAVS3/); DIMENSIONS_MOD gives the array size JSON.

The method

Inside the grid (IBEGIN < $i$ < IEND): for each model point, find the four grid points around its radius, get Lagrange weights $b_k$, and $$\Delta X_i = \sum_{k=1}^{4} b_k\,\Delta X(r_k), \qquad X_i \leftarrow \min(X_i + \Delta X_i,\; 1 - Z_i - X_{3,i}),$$ where $X_3$ is $^3$He (HCOMP(4,·)). With metal diffusion, $Z$ changes the same way, and every metal species (rows 5–11) is scaled by $Z_\mathrm{new}/Z_\mathrm{old}$. Helium takes the rest: $Y = 1 - X - Z - X_3$. The scaling divides by the old $Z$ (suspected bug 1). The line IF (J .EQ. 0) PRINT*, 'line 47 etm' (line 69 ↗) cannot fire, since J starts at JMIN ≥ 2 (unused?); nor can IF (K0 .EQ. 0) K0=1 (lines 76–77 ↗), since NPT is at least 4 (unused?).
Below the grid ($i \le$ IBEGIN) the change at the first grid point is applied; above it ($i \ge$ IEND), the change at the last. The points IBEGIN and IEND are themselves the grid's end points: MODEL_TO_EQUAL puts the first grid point at HRU(IBEGIN) and spaces the grid so that the last falls at HRU(IEND), with the model's values there (settled from the code, Nov 9, 2026: lines 42–44 ↗, lines 113–116 ↗; model_to_equal.f, line 73 ↗, model_to_equal.f, lines 195–198 ↗, model_to_equal.f, lines 205–206 ↗, model_to_equal.f, lines 266–267 ↗). The comment explains the guards: "because metals sink, and hydrogen rises, the failsafes are opposite" ($X \ge 0$ and $Z \le Z_\mathrm{max}$ below; $X \le X_\mathrm{max}$ and $Z \ge 0$ above).
Units: $r \leftarrow r/c_r$, $T \leftarrow T/c_T$, $m \leftarrow m/c_m$, HQPR $\leftarrow$ HQPR$\cdot c_r$, $\Delta t \leftarrow \Delta t\,c_t$, $M \leftarrow M/c_m$.

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) Division by $Z$.

(line 52 ↗; line 103 ↗; line 122 ↗) With metal diffusion, each shell's metal species are scaled by ZZ2 = ZZ/HCOMP(3,I), the new $Z$ over the old; a shell with $Z = 0$ divides by zero, and rows 5–11 become NaN (or infinite). What follows: A metal-free model (Population III) cannot be run with LDIFZ. Only with LDIFZ (off by default) and LNEWDIF false, in a shell with $Z = 0$.

A possible fix (the same at lines 103 and 122):

equal_to_model.f

@@ lines 49–55 @@
          DO I = IBEGIN,1,-1
             ZMAX = 1.0D0 - HCOMP(1,I) - HCOMP(4,I)
             ZZ=MIN(HCOMP(3,I)+EZ(1),ZMAX)
-            ZZ2 = ZZ/HCOMP(3,I)
+C 11/09/26 KC via Claude
+C code changed to leave the metal species unscaled where Z was zero.
+            IF(HCOMP(3,I).GT.0.0D0)THEN
+               ZZ2 = ZZ/HCOMP(3,I)
+            ELSE
+               ZZ2 = 1.0D0
+            ENDIF
             HCOMP(3,I) = ZZ
             DO J = 5,11
                HCOMP(J,I) = ZZ2*HCOMP(J,I)
@@ lines 100–106 @@
             DZMOD = FACINTERP(1)*EZ(K0)+FACINTERP(2)*EZ(K0+1)+
      *             FACINTERP(3)*EZ(K0+2)+FACINTERP(4)*EZ(K0+3)
             ZZ = MIN(HCOMP(3,I)+DZMOD,ZMAX)
-            ZZ2 = ZZ/HCOMP(3,I)
+C 11/09/26 KC via Claude
+C code changed to leave the metal species unscaled where Z was zero.
+            IF(HCOMP(3,I).GT.0.0D0)THEN
+               ZZ2 = ZZ/HCOMP(3,I)
+            ELSE
+               ZZ2 = 1.0D0
+            ENDIF
             HCOMP(3,I)=ZZ
             DO J = 5,11
                HCOMP(J,I) = ZZ2*HCOMP(J,I)
@@ lines 119–125 @@
          ZMIN = 0.0D0
          DO I = IEND,M
             ZZ = MAX(HCOMP(3,I)+EZ(NPT),ZMIN)
-            ZZ2 = ZZ/HCOMP(3,I)
+C 11/09/26 KC via Claude
+C code changed to leave the metal species unscaled where Z was zero.
+            IF(HCOMP(3,I).GT.0.0D0)THEN
+               ZZ2 = ZZ/HCOMP(3,I)
+            ELSE
+               ZZ2 = 1.0D0
+            ENDIF
             HCOMP(3,I) = ZZ
             DO J = 5,11
                HCOMP(J,I) = ZZ2*HCOMP(J,I)

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 ↗), and in addition a shell with $Z = 0$; no case has a metal-free model. Intent: none found.

Definitions

subroutineEQUAL_TO_MODEL

lines 2–144 · GitHub

Interpolates the settling changes in $X$ (and $Z$) from the equally spaced grid to the model points with four-point Lagrange weights, enforces the composition limits, and converts the model back to physical units.

Uses modules CONFAC_MOD, DIMENSIONS_MOD, GRAVEZ_MOD, NMLPAR_MOD

Called by GRSETT

Flowchart of its calls, in order (1)

for I=IBEGIN+1,IEND-1 (DO 20), lines 65–112

INTRP2 91

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE.

1C $$$$$$ 2 SUBROUTINE EQUAL_TO_MODEL(DT,ER,EX,IBEGIN,IEND,NPT,HCOMP,HQPR, 3 * * HRU,HS1,HTU,M,STOT) 4C 10/15/26 MHP via Claude 5C code changed to take JSON from DIMENSIONS_MOD instead of declaring 6C PARAMETER(JSON=5000) here. 7 USE DIMENSIONS_MOD 8C 10/14/26 MHP via Claude 9C code changed to convert COMMON/GRAVEZ/ into GRAVEZ_MOD. 10 USE GRAVEZ_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/CONFAC/ into CONFAC_MOD. 13 USE CONFAC_MOD 14C 10/14/26 MHP via Claude 15C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 16 USE NMLPAR_MOD 17C 10/20/26 MHP via Claude 18C code changed to replace implicit typing with IMPLICIT NONE and 19C explicit declarations. 20 IMPLICIT NONE 21 REAL*8 DT,DXMOD,DZMOD,ER,EX,FACINTERP,HCOMP,HQPR,HRU,HS1,HTU, 22 * * RADMOD,STOT,TABLER,XMAX,XMINIM,ZMAX,ZMIN,ZZ,ZZ2 23 INTEGER I,IBEGIN,IEND,J,JMIN,K,K0,M,NPT 24C 10/12/26 MHP via Claude 25C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 26C namelist variables that the physics routines read from 27 COMMON 27C blocks which also held derived or working members. Those blocks keep 28C only their non-namelist members; five of them -- NEWMX, VMULT, 29C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 30C Derivation-only namelist variables from the same blocks went to 31C COMMON/NMLBYPASS/ instead, not here. 32C 33C Declared identically in all 121 files that hold any of these members. 34C Generated from one member list, so every copy is byte-identical; the 35C linker's 'size of symbol changed' warning is the check that they are. 36 DIMENSION ER(JSON),EX(JSON),HCOMP(15,JSON),HQPR(JSON),HRU(JSON), 37 * * HS1(JSON),HTU(JSON),TABLER(4),FACINTERP(4) 38 SAVE 39C TRANSFORM BACK TO ORIGINAL GRID OF MODEL POINTS FROM EQUALLY 40C SPACED GRID. 41 XMINIM = 0.0D0 42 DO 10 I = IBEGIN,1,-1 43 HCOMP(1,I)=MAX(HCOMP(1,I) + EX(1),XMINIM) 44 10 10 CONTINUE 45C MHP 3/94 ADDED METAL DIFFUSION 46C NOTE THAT BECAUSE METALS SINK, AND HYDROGEN RISES, THE FAILSAFES 47C ARE OPPOSITE (GUARDING AGAINST NEGATIVE X AND Z>1 RESPECTIVELY). 48 IF(LDIFZ)THEN 49 DO I = IBEGIN,1,-1 50 ZMAX = 1.0D0 - HCOMP(1,I) - HCOMP(4,I) 51 ZZ=MIN(HCOMP(3,I)+EZ(1),ZMAX) 52 ZZ2 = ZZ/HCOMP(3,I) 53 HCOMP(3,I) = ZZ 54 DO J = 5,11 55 HCOMP(J,I) = ZZ2*HCOMP(J,I) 56 END DO 57 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 58 END DO 59 ELSE 60 DO I = IBEGIN,1,-1 61 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 62 END DO 63 ENDIF 64 JMIN=2 65 DO 20 I=IBEGIN+1,IEND-1 66 DO 30 J=JMIN,NPT 67 68C FIND EQUALLY SPACED GRID POINTS CLOSEST TO THE MODEL POINT. 69 IF (J .EQ. 0) PRINT*, 'line 47 etm' 70 IF(ER(J).GE.HRU(I))THEN 71 72C ENSURE THAT FIRST INTERP. POINT NO LESS THAN FIRST EQUALLY SPACED POINT. 73 K0 = MAX(J-2,1) 74C ENSURE THAT LAST INTERP. POINT NO GREATER THAN LAST EQUALLY SPACED POINT. 75 K0 = MIN(K0,NPT-3) 76C JVS fix for NPT = 3? 77 IF (K0 .EQ. 0) K0=1 78 JMIN=J 79 GOTO 40 80 ENDIF 81 30 30 CONTINUE 82 K0 = NPT-3 83 JMIN=NPT 84 40 40 CONTINUE 85 86 DO 50 K=1,4 87 TABLER(K)=ER(K0+K-1) 88 50 50 CONTINUE 89 RADMOD=HRU(I) 90C FIND 4 POINT LAGRANGIAN INTERPOLATION FACTORS. 91 CALL INTRP2(TABLER,FACINTERP,RADMOD) 92C PERFORM 4 POINT LAGRANGIAN INTERPOLATION FOR CHANGE IN X. 93 DXMOD = FACINTERP(1)*EX(K0)+FACINTERP(2)*EX(K0+1)+ 94 * * FACINTERP(3)*EX(K0+2)+FACINTERP(4)*EX(K0+3) 95 XMAX = 1.0D0 - HCOMP(3,I) - HCOMP(4,I) 96 HCOMP(1,I)=MIN(HCOMP(1,I) + DXMOD,XMAX) 97C MHP 3/94 ADDED METAL DIFFUSION 98 IF(LDIFZ)THEN 99 ZMAX = 1.0D0 - HCOMP(1,I) - HCOMP(4,I) 100 DZMOD = FACINTERP(1)*EZ(K0)+FACINTERP(2)*EZ(K0+1)+ 101 * * FACINTERP(3)*EZ(K0+2)+FACINTERP(4)*EZ(K0+3) 102 ZZ = MIN(HCOMP(3,I)+DZMOD,ZMAX) 103 ZZ2 = ZZ/HCOMP(3,I) 104 HCOMP(3,I)=ZZ 105 DO J = 5,11 106 HCOMP(J,I) = ZZ2*HCOMP(J,I) 107 END DO 108 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 109 ELSE 110 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 111 ENDIF 112 20 20 CONTINUE 113 DO 60 I = IEND,M 114 XMAX = 1.0D0 - HCOMP(3,I) - HCOMP(4,I) 115 HCOMP(1,I)=MIN(HCOMP(1,I) + EX(NPT),XMAX) 116 60 60 CONTINUE 117C MHP 3/94 ADDED METAL DIFFUSION 118 IF(LDIFZ)THEN 119 ZMIN = 0.0D0 120 DO I = IEND,M 121 ZZ = MAX(HCOMP(3,I)+EZ(NPT),ZMIN) 122 ZZ2 = ZZ/HCOMP(3,I) 123 HCOMP(3,I) = ZZ 124 DO J = 5,11 125 HCOMP(J,I) = ZZ2*HCOMP(J,I) 126 END DO 127 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 128 END DO 129 ELSE 130 DO I = IEND,M 131 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 132 END DO 133 ENDIF 134 DO 70 I=1,M 135 136 HRU(I)=HRU(I)/CON_RAD 137 HTU(I)=HTU(I)/CON_TEMP 138 HS1(I)=HS1(I)/CON_MASS 139 HQPR(I)=HQPR(I)*CON_RAD 140 70 70 CONTINUE 141 DT=DT*CON_TIME 142 STOT=STOT/CON_MASS 143 RETURN 144 END

1"""Python version of equal_to_model.f (YREC 7be4133): the settling 2changes from the equally spaced grid back to the model points.""" 3 4import numpy as np 5 6from pyrec import confac_mod 7from pyrec import gravez_mod 8from pyrec import nmlpar_mod 9from pyrec.intrp2 import intrp2 10 11 12def _scale_metals(hcomp: np.ndarray, i: int, zz: float) -> None: 13 """Set Z of shell i to zz, scale the metal species (rows 5-11) by 14 Z_new/Z_old, and give helium the rest (Fortran lines 52–57, 103–108, 15 122–127). 16 17 """ 18 19 # Suspected bug (see the page): a shell with Z = 0 divides by zero. 20 zz2 = zz / hcomp[3, i] 21 hcomp[3, i] = zz 22 for j in range(5, 11 + 1): 23 hcomp[j, i] = zz2 * hcomp[j, i] 24 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 25 26 27def equal_to_model(dt: float, er: np.ndarray, ex: np.ndarray, ibegin: int, iend: int, npt: int, hcomp: np.ndarray, 28 hqpr: np.ndarray, hru: np.ndarray, hs1: np.ndarray, htu: np.ndarray, m: int, 29 stot: float) -> tuple[float, float]: 30 """Apply the changes in X (and, with nmlpar_mod.ldifz, in Z) found on 31 the equally spaced grid to the model points, by four-point Lagrange 32 interpolation inside the grid and the end values outside it, within 33 the composition limits; then convert the model back from the 34 solver's units. 35 36 Translates: equal_to_model.f#EQUAL_TO_MODEL (7be4133) 37 38 Arguments: er, ex (grid radii and changes in X, npt points); ibegin, 39 iend (the model points bounding the grid); m (number of shells). 40 Arrays changed in place: hcomp, hqpr, hru, hs1, htu. 41 42 Returns (dt, stot), back in physical units. 43 44 """ 45 46 tabler = np.zeros(4 + 1) 47 facinterp = np.zeros(4 + 1) 48 ez = gravez_mod.ez 49 50 # Below the grid: the change at the first grid point (Fortran lines 51 # 41–44). 52 xminim = 0.0 53 for i in range(ibegin, 0, -1): # label 10 54 hcomp[1, i] = max(hcomp[1, i] + ex[1], xminim) 55 56 # Metals sink and hydrogen rises, so the failsafes are opposite 57 # (Fortran lines 48–63). 58 if nmlpar_mod.ldifz: 59 for i in range(ibegin, 0, -1): 60 zmax = 1.0 - hcomp[1, i] - hcomp[4, i] 61 zz = min(hcomp[3, i] + ez[1], zmax) 62 _scale_metals(hcomp, i, zz) 63 else: 64 for i in range(ibegin, 0, -1): 65 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 66 67 # Inside the grid: four-point Lagrange interpolation (Fortran lines 68 # 64–112). 69 jmin = 2 70 for i in range(ibegin + 1, iend): # label 20 71 for j in range(jmin, npt + 1): # label 30 72 if j == 0: # never true 73 print(' line 47 etm') 74 if er[j] >= hru[i]: 75 k0 = max(j - 2, 1) 76 k0 = min(k0, npt - 3) 77 if k0 == 0: 78 k0 = 1 79 jmin = j 80 break # GOTO 40 81 else: 82 k0 = npt - 3 83 jmin = npt 84 # label 40 85 for k in range(1, 4 + 1): # label 50 86 tabler[k] = er[k0 + k - 1] 87 radmod = hru[i] 88 intrp2(tabler, facinterp, radmod) 89 dxmod = (facinterp[1] * ex[k0] + facinterp[2] * ex[k0 + 1] 90 + facinterp[3] * ex[k0 + 2] + facinterp[4] * ex[k0 + 3]) 91 xmax = 1.0 - hcomp[3, i] - hcomp[4, i] 92 hcomp[1, i] = min(hcomp[1, i] + dxmod, xmax) 93 if nmlpar_mod.ldifz: 94 zmax = 1.0 - hcomp[1, i] - hcomp[4, i] 95 dzmod = (facinterp[1] * ez[k0] + facinterp[2] * ez[k0 + 1] 96 + facinterp[3] * ez[k0 + 2] 97 + facinterp[4] * ez[k0 + 3]) 98 zz = min(hcomp[3, i] + dzmod, zmax) 99 _scale_metals(hcomp, i, zz) 100 else: 101 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 102 103 # Above the grid: the change at the last grid point (Fortran lines 104 # 113–116). 105 for i in range(iend, m + 1): # label 60 106 xmax = 1.0 - hcomp[3, i] - hcomp[4, i] 107 hcomp[1, i] = min(hcomp[1, i] + ex[npt], xmax) 108 109 # Fortran lines 118–133. 110 if nmlpar_mod.ldifz: 111 zmin = 0.0 112 for i in range(iend, m + 1): 113 zz = max(hcomp[3, i] + ez[npt], zmin) 114 _scale_metals(hcomp, i, zz) 115 else: 116 for i in range(iend, m + 1): 117 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 118 119 # Back to physical units (Fortran lines 134–142). 120 for i in range(1, m + 1): # label 70 121 hru[i] = hru[i] / confac_mod.con_rad 122 htu[i] = htu[i] / confac_mod.con_temp 123 hs1[i] = hs1[i] / confac_mod.con_mass 124 hqpr[i] = hqpr[i] * confac_mod.con_rad 125 dt = dt * confac_mod.con_time 126 stot = stot / confac_mod.con_mass 127 128 return dt, stot