rotmix.f

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

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

Variables on this page (14)

cgl, cln, csecyr, delm, dt_gs, ishort, json, ldify, lexcom, lnewdif, sdel, tcut, xmin, ymin

Purpose

Nuclear burning and gravitational settling for one step of the rotational-mixing loop. Its header comment lists the steps: (1) use KEMCOM to solve implicitly for the abundances in radiative regions, shell by shell; (2) solve for the burning in each convection zone as a whole. It then applies the settling of helium and renormalizes the composition.

How it fits

Who calls it: MIDMOD (with the temperatures at the current time, HTM, "for thermodynamic consistency," MHP 6/00) and BURSMIX, in their loops of substeps.
Inputs: the step DELTS (s), the composition HCOMP(15, ·), shell masses HS2, temperatures HT ($\log T$), densities HD (also used for the rates), and the zone lists MRZONE (radiative) and MXZONE (convective); for the settling, HS, HR, HP, HS1, HSTOT and the convection flags LC. The zone lists come from CONVEC, whose end of the surface radiative zone is misplaced (convec.f, suspected bug 1).
The rates: computed here by RATES on the model passed in (HD, HT, and the substep's own composition), as MIX does: the rates HR1…HR13 and the branchings HF1, HF2, set to zero in shells with $\log T \le$ TCUT(1).
What it calls: RATES, KEMCOM, then either MICRODIFF (when LNEWDIF is true) or GRSETT (the older Bahcall–Loeb settling).
Output: HCOMP, updated in place.
Modules, other: NMLPAR_MOD (the namelist parameters LDIFY, which this routine may switch off, DT_GS, XMIN, YMIN, LNEWDIF, TCUT(1) and LEXCOM: 15 species or 11), MDPHY_MOD (DELM), SCRTCH_MOD (SDEL, temporarily overwritten), CONSTANTS_MOD, IOUNITS_MOD (ISHORT), DIMENSIONS_MOD (JSON).

The method

Burning. For each radiative shell above the cut $\log T >$ TCUT(1) (the loop stops at the first shell below it, assuming $T$ falls outward), KEMCOM is called for the shell alone; for each convection zone, for the whole zone at once, which also homogenizes it. The step is passed in years, DDAGE = DELTS/CSECYR.

The burning loop's early exit skips nothing hot. GOTO 45 leaves both the shell loop and the zone loop at the first shell cooler than TCUT(1) (lines 118–130 ↗), so a radiative zone listed after it is skipped; but CONVEC lists the radiative zones from the center out (convec.f, lines 185–205 ↗, with the slip of CONVEC's suspected bug 1 in the last zone's upper edge), and the temperature falls outward, so every later zone is cooler too. MIX exits its loop in the same way (mix.f, line 252 ↗).
Settling (if LDIFY). It is skipped once the central hydrogen is below XMIN (and LDIFY is then switched off for good), if the star is fully convective, or if $Y <$ YMIN everywhere. It stops the run (exit status 1) if there is no surface convection zone. The outer point JMAX is the top of the last radiative zone; the comment says it is moved down to the first shell with $Y >$ YMIN, but the search moves only J, so JMAX stays the top of the last radiative zone (suspected bug 1). The settling time there is taken from Michaud et al. (1984): $$t_\mathrm{settle} = 4.348\times10^{21}\;\frac{f_M}{T^{3/2}}\ \text{yr},\qquad f_M = \frac{M_* - m(\texttt{JMAX})}{M_*},$$ with $T$ in K at JMAX (the code multiplies by CSECYR to get TSCALE in seconds, the unit of DELTS), and the step is divided into NSTEP substeps of at most DT_GS $\times\,t_\mathrm{settle}$ (one more than INT(DELTS/DTMAX) unless DELTS is an exact multiple of DTMAX). The pressure gradient passed to the settling routine is $$\frac{d\ln P}{dr} = -\frac{G m \rho}{P r^2}\quad(\texttt{HQPR}).$$ During the settling, SDEL(2,·) holds DELM and is restored afterwards.
Renormalization. Every abundance is clipped to $[0, 1]$; $Z$ (HCOMP(3)) is capped at $1 - X - X_3$, and $Y$ is set to $1 - X - Z - X_3$, where $X_3$ is ³He, the fourth species HCOMP(4,·) (as STARIN fills it).

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. (never occurred) The outer point for the settling step is not moved down.

(lines 156–160 ↗; lines 174–180 ↗; lines 191–194 ↗) The comment says the outer boundary is the surface convection zone or "the first point where the helium abundance rises above" YMIN. Loop 140 finds that point as J, but JMAX is never set to it, so FM and TSCALE, and with them the number of settling substeps, are always taken at the top of the last radiative zone. Where $Y \le$ YMIN there, the step limit is set at a point the comment means to skip. Only with LDIFY (off by default); the settling itself does not use JMAX.

A possible fix:

rotmix.f

@@ lines 178–183 @@
 C   Y<YMIN FOR THE WHOLE STAR IF THE CODE GETS HERE.
          GOTO 170
   150    CONTINUE
+C 11/09/26 KC via Claude
+C code changed to move the outer point down to the first shell with
+C Y > YMIN, as the comment above says.
+         JMAX = J
          STOT=EXP(CLN*HSTOT)
          DO 130 I = 1,M
             TEMP(I) = SDEL(2,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): The outer point was searched in the two rotating test cases with settling, 2,217 times in testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS and 1,658 in testsuite/Test_solarGS98_difrotmixper18p4d1p0, at shells 8 to 527, and the first shell with $Y >$ YMIN was always the top of the last radiative zone itself, so J equaled JMAX and the missing assignment changed nothing. It matters only where the helium abundance there has fallen to YMIN ($10^{-3}$ in every namelist), which no case reaches. Intent: MIX, the same code for runs without rotation, uses J (mix.f, line 441 ↗, mix.f, line 444 ↗), so this looks like a slip.

2. (not reached) The stop message names the wrong routine.

(lines 163–169 ↗) The stop for no surface convection zone prints "STOPPED IN SUBROUTINE MIX" (copied from MIX), though it stops in ROTMIX. Only with LDIFY (off by default).

A possible fix:

rotmix.f

@@ lines 163–170 @@
          IF(.NOT.LC(M))THEN
             WRITE(*,911)
             WRITE(ISHORT,911)
+C 11/09/26 KC via Claude
+C code changed to name ROTMIX in the stop message.
   911       FORMAT(1X,'NO SURFACE CZ - DIFFUSION NOT MEANINGFUL'/
-     *             'STOPPED IN SUBROUTINE MIX')
+     *             'STOPPED IN SUBROUTINE ROTMIX')
 C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38).
             STOP 1
          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. No case run stopped here: The rotating cases with settling (testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, testsuite/Test_solarGS98_difrotmixper18p4d1p0) have a surface convection zone, and the 2 and 3 $M_\odot$ rotating cases have LDIFY off. A rotating run with LDIFY for a star without a surface convection zone would reach it, and stop. Intent: The text was copied from MIX, which since 1992 only warns and goes on (mix.f, lines 428–433 ↗: "DBG 2/92 CHANGED STOP TO JUST A WARNING MESSAGE, EXECUTION CONTINUES").

Definitions

subroutineROTMIX

lines 5–251 · GitHub

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

Called by BURSMIX, MIDMOD

Flowchart of its calls, in order (5)

for II = 1,M, lines 87–110

else, lines 104–109

RATES 105

for K = 1,NRZONE (DO 40), lines 118–129

for J = MRZONE(K,1),MRZONE(K,2) (DO 30), lines 119–128

KEMCOM 124

for K = 1,NZONE (DO 50), lines 135–150

KEMCOM 138

if (LDIFY), lines 151–237

for NS = 1,NSTEP (DO 160), lines 224–232

if (LNEWDIF), lines 227–228

MICRODIFF 228

else, lines 229–231

GRSETT 230

Compiler warnings (1)

gfortran -Wall -Wextra, compiling only; nothing was run.

line 207: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. The two debug writes to unit 69 are commented out (settles former question 1); COMMON blocks to modules (the unused /CTLIM/ dropped), IMPLICIT NONE.

Changed since a8a5950. The rates are computed by RATES on the model passed in, instead of copied from HCOMPM in BURN_MOD (removed since a8a5950), which held the rates MIX found on the previous Newton iterate's end-of-step structure, so that the midpoint substep from MIDMOD burned with end-of-step rates (this settles former question 3, since the copies are gone); the substep count tests MOD(DELTS,DTMAX) (this settles former question 1); the stop for no surface convection zone exits with status 1, the code for run-time errors (KNOWN_ISSUES.md, issue 38). The later questions are renumbered.

1C 1) USE KEMCOM TO IMPLICITLY SOLVE FOR ABUNDANCES IN RADIATIVE REGIONS. 2C 2) SOLVE FOR BURNING IN EACH CONVECTION ZONE. 3 4C SUBROUTINE ROTMIX(DELTS,HCOMP,HS2,HT,ITLVL,M,MRZONE,MXZONE, ! KC 2025-05-31 5 SUBROUTINE ROTMIX(DELTS,HCOMP,HS2,HT,M,MRZONE,MXZONE, 6 * * NRZONE,NZONE 7 * * ,HSTOT,HD,HS,HR,HP,LC,HS1) 8C 10/15/26 MHP via Claude 9C code changed to take JSON from DIMENSIONS_MOD instead of declaring 10C PARAMETER(JSON=5000) here. 11 USE DIMENSIONS_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/SCRTCH/ into SCRTCH_MOD. 14 USE SCRTCH_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/MDPHY/ into MDPHY_MOD. 17 USE MDPHY_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 20 USE NMLPAR_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 23 USE CONSTANTS_MOD 24C 10/14/26 MHP via Claude 25C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 26 USE IOUNITS_MOD 27C 10/20/26 MHP via Claude 28C code changed to replace implicit typing with IMPLICIT NONE and 29C explicit declarations. 30 IMPLICIT NONE 31 REAL*8 DDAGE,DELTS,DT,DTMAX,FM,HCOMP,HD,HF1,HF2,HP,HQPR,HR,HR1, 32 * * HR10,HR11,HR12,HR13,HR2,HR3,HR4,HR5,HR6,HR7,HR8,HR9,HS,HS1, 33 * * HS2,HSTOT,HT,STOT,TEMP,TSCALE 34 INTEGER I,IBEGIN,IEND,II,J,JMAX,K,M,MRZONE,MXZONE,NRZONE,NS,NSPEC, 35 * * NSTEP,NZONE 36 LOGICAL*4 LC,LCZ 37 38C 9/25/26 MHP via Claude 39C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 40C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 41C 9/24/26 MHP via Claude 42C code changed to remove dead CMIXL from COMMON/CONST3/; it is now 43C passed as an explicit argument where actually used 44C 10/12/26 MHP via Claude 45C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 46C namelist variables that the physics routines read from 27 COMMON 47C blocks which also held derived or working members. Those blocks keep 48C only their non-namelist members; five of them -- NEWMX, VMULT, 49C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 50C Derivation-only namelist variables from the same blocks went to 51C COMMON/NMLBYPASS/ instead, not here. 52C 53C Declared identically in all 121 files that hold any of these members. 54C Generated from one member list, so every copy is byte-identical; the 55C linker's 'size of symbol changed' warning is the check that they are. 56C 10/13/26 MHP via Claude 57C code changed to drop COMMON/CTLIM/, declared here but never used. 58c GES 6/15 INCLUDED COMMON BLOCK FOR NEW DIFFUSION ROUTINES. 59C 9/25/26 MHP via Claude 60C code changed to consolidate all I/O logical unit numbers, previously 61C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 62C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 63C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 64C COMMON/IOUNITS/ 65 DIMENSION HS2(JSON),HT(JSON),HCOMP(15,JSON) 66 DIMENSION HR1(JSON),HR2(JSON),HR3(JSON),HR4(JSON), 67 * * HR5(JSON),HR6(JSON),HR7(JSON),HR8(JSON),HR9(JSON),HR10(JSON), 68 * * HR11(JSON),HR12(JSON),HR13(JSON),HF1(JSON),HF2(JSON) 69 DIMENSION MRZONE(13,2),MXZONE(12,2), 70 * * HD(JSON),HS(JSON),HR(JSON),HP(JSON),HS1(JSON),HQPR(JSON), 71 * * LC(JSON),TEMP(JSON),LCZ(JSON) 72 SAVE 73C NSPEC IS THE NUMBER OF SPECIES BEING TRACKED. 74 IF(LEXCOM)THEN 75 NSPEC = 15 76 ELSE 77 NSPEC = 11 78 ENDIF 79C DDAGE IS THE TIMESTEP IN YEARS. 80 DDAGE=DELTS/CSECYR 81C 10/29/26 MHP via Claude 82C code changed to compute the burning rates on the model passed in (HD, 83C HT and the sub-step's own composition), as MIX does, instead of copying 84C HCOMPM. HCOMPM held the rates MIX found on the previous Newton 85C iterate's end-of-step structure, so the midpoint sub-step from MIDMOD 86C burned with end-of-step rates, contrary to the MHP 6/00 intent there. 87 DO II = 1,M 88 IF(HT(II).LE.TCUT(1))THEN 89 HR1(II) = 0.0D0 90 HR2(II) = 0.0D0 91 HR3(II) = 0.0D0 92 HR4(II) = 0.0D0 93 HR5(II) = 0.0D0 94 HR6(II) = 0.0D0 95 HR7(II) = 0.0D0 96 HR8(II) = 0.0D0 97 HR9(II) = 0.0D0 98 HR10(II) = 0.0D0 99 HR11(II) = 0.0D0 100 HR12(II) = 0.0D0 101 HR13(II) = 0.0D0 102 HF1(II) = 0.0D0 103 HF2(II) = 0.0D0 104 ELSE 105 CALL RATES(HD(II),HT(II),HCOMP(1,II),HCOMP(2,II), 106 * * HCOMP(4,II),HCOMP(5,II),HCOMP(6,II),HCOMP(7,II), 107 * * HCOMP(9,II),HCOMP(11,II),II,HR1,HR2,HR3,HR4,HR5,HR6, 108 * * HR7,HR8,HR9,HR10,HR11,HR12,HR13,HF1,HF2) 109 ENDIF 110 END DO 111C 112C NOW IMPLICITLY SOLVE FOR THE NEW ABUNDANCES AT THE END OF THE 113C TIMESTEP. THIS IS DONE SHELL BY SHELL FOR RADIATIVE REGIONS, 114C AND FOR EACH CONVECTION ZONE AS A UNIT. 115C 116C RADIATIVE ZONES. 117C 118 DO 40 K = 1,NRZONE 119 DO 30 J = MRZONE(K,1),MRZONE(K,2) 120C EXIT LOOP ONCE T DROPS BELOW NUCLEAR REACTION T CUTOFF 121 IF(HT(J).LE.TCUT(1)) GOTO 45 122 IBEGIN = J 123 IEND = J 124 CALL KEMCOM(HT,IBEGIN,IEND,HR1,HR2,HR3,HR4,HR5,HR6,HR7, 125C * HR8,HR9,HR10,HR11,HR12,HR13,HF1,HS2,HCOMP, 126C * DDAGE,ITLVL) ! KC 2025-05-31 127 * * HR8,HR10,HR11,HR12,HF1,HS2,HCOMP,DDAGE) 128 30 30 CONTINUE 129 40 40 CONTINUE 130 45 45 CONTINUE 131C 132C CONVECTION ZONES. 133C NOTE KEMCOM ALSO AUTOMATICALLY HOMOGENIZE CONVECTION ZONES. 134C 135 DO 50 K = 1,NZONE 136 IBEGIN = MXZONE(K,1) 137 IEND = MXZONE(K,2) 138 CALL KEMCOM(HT,IBEGIN,IEND,HR1,HR2,HR3,HR4,HR5,HR6,HR7, 139C * HR8,HR9,HR10,HR11,HR12,HR13,HF1,HS2,HCOMP, 140C * DDAGE,ITLVL) ! KC 2025-05-31 141 * * HR8,HR10,HR11,HR12,HF1,HS2,HCOMP,DDAGE) 142 50 50 CONTINUE 143C 144C MICROSCOPIC DIFFUSION OF HELIUM. 145C 146C***BC 5/92 ROTMIX MODIFIED TO INCLUDE CALL TO GRAVITATIONAL SETTLING 147C ROUTINE USING THE BAHCALL AND LOEB METHOD. 148C FIRST DEFINE VARIABLES NEEDED FOR SETTLING - 149C HQPR=VECTOR OF D LN P/DR. 150C STOT=TOTAL STELLAR MASS(UNLOGGED). 151 IF(LDIFY)THEN 152 IF(HCOMP(1,1).LT.XMIN)THEN 153 LDIFY=.FALSE. 154 GOTO 170 155 ENDIF 156C MHP 6/90 CHANGE ADDED : THE TIMESTEP FOR SETTLING IS RESTRICTED TO 157C A FRACTION OF THE TIMESCALE FOR SETTLING AT THE OUTER BOUNDARY. 158C THE OUTER BOUNDARY IS EITHER THE SURFACE CONVECTION ZONE OR THE 159C FIRST POINT WHERE THE HELIUM ABUNDANCE RISES ABOVE A USER-SPECIFIED 160C MINIMUM VALUE (YMIN). 161C 162C LOCATE OUTER BOUNDARY. 163 IF(.NOT.LC(M))THEN 164 WRITE(*,911) 165 WRITE(ISHORT,911) 166 911 911 FORMAT(1X,'NO SURFACE CZ - DIFFUSION NOT MEANINGFUL'/ 167 * * 'STOPPED IN SUBROUTINE MIX') 168C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 169 STOP 1 170 ENDIF 171C 7/92 MHP STATEMENT ADDED FOR EXIT IF OVERSHOOT CAUSES A FULLY CONVECTIVE 172C CASE. AVOIDS OCCASIONAL PRE-MS CRASH IN FIRST RADIATIVE MODEL. 173 IF(MXZONE(NZONE,1).LE.2 .AND. MXZONE(NZONE,2).EQ.M)GOTO 170 174 JMAX = MRZONE(NRZONE,2) 175 DO 140 J = JMAX,1,-1 176 IF(HCOMP(2,J).GT.YMIN)GOTO 150 177 140 140 CONTINUE 178C Y<YMIN FOR THE WHOLE STAR IF THE CODE GETS HERE. 179 GOTO 170 180 150 150 CONTINUE 181 STOT=EXP(CLN*HSTOT) 182 DO 130 I = 1,M 183 TEMP(I) = SDEL(2,I) 184 SDEL(2,I) = DELM(I) 185 HQPR(I)=-EXP(CLN*(HD(I)+CGL+HS(I)-2.0D0*HR(I)-HP(I))) 186 130 130 CONTINUE 187C ***BC 6/92 only check for timestep cutting if JMAX is large 188C 189 IF(JMAX.GE.2) THEN 190C FM IS THE MASS FRACTION ABOVE THE OUTER POINT. 191 FM = (STOT-HS1(JMAX))/STOT 192C TSCALE IS THE TIMESCALE FOR SETTLING OF HELIUM AT THE OUTER 193C BOUNDARY (MICHAUD ET AL 1984, APJ V.282,P.206) 194 TSCALE = 4.348D21*CSECYR*FM/EXP(CLN*1.5D0*HT(JMAX)) 195C RESTRICT TIMESTEP TO THE MINIMUM OF THE MODEL TIMESTEP AND 196C A USER SPECIFIED FRACTION (DT_GS) OF THE SETTLING TIMESCALE. 197C 9/25/26 MHP via Claude 198C code changed to comment out debug write to unit 69, which was never 199C opened anywhere (silently produced a stray fort.69 file). Left in 200C place, commented, rather than deleted. 201C write(69,*) 'JMAX=',JMAX,' FM= ',fm,' TSCALE=',TSCALE 202 DTMAX = DT_GS*TSCALE 203 NSTEP = INT(DELTS/DTMAX) 204C 10/30/26 MHP via Claude 205C code changed to MOD(DELTS,DTMAX); the arguments were reversed, which 206C took one sub-step too many when DELTS was an exact multiple of DTMAX. 207 IF(MOD(DELTS,DTMAX).NE.0.0D0.OR.NSTEP.EQ.0)NSTEP=NSTEP+1 208C write(69,*)'DTMAX=',dtmax,' DELTS=',delts,' NSTEP=',nstep 209 DT = DELTS/DFLOAT(NSTEP) 210 ELSE 211 NSTEP=1 212 DT = DELTS 213 ENDIF 214 DO J = 1,NRZONE 215 DO I = MRZONE(J,1),MRZONE(J,2) 216 LCZ(I) = .FALSE. 217 END DO 218 END DO 219 DO J = 1,NZONE 220 DO I = MXZONE(J,1),MXZONE(J,2) 221 LCZ(I) = .TRUE. 222 END DO 223 END DO 224 DO 160 NS = 1,NSTEP 225C PERFORM GRAVITATIONAL SETTLING. IF LNEWDIF = TRUE, USE THE NEW ROUTINES 226C IN MICRODIFF. ELSE, USE THE OLD ROUTINES IN GRSETT. 227 IF(LNEWDIF)THEN 228 CALL MICRODIFF(DT,HCOMP,HQPR,HR,HD,HS1,HT,LCZ,M,STOT) 229 ELSE 230 CALL GRSETT(DT,HCOMP,HQPR,HR,HD,HS1,HT,LCZ,M,STOT) 231 ENDIF 232 160 160 CONTINUE 233 DO I = 1,M 234 SDEL(2,I) = TEMP(I) 235 END DO 236 170 170 CONTINUE 237 ENDIF 238C END OF GRAVITATIONAL SETTLING 239C 240C RENORMALIZE COMPOSITION TO GUARD AGAINST ANOMALIES (I.E. SMALL NEGATIVE 241C ABUNDANCES...). 242 DO 180 I = 1,M 243 DO 175 J = 1,NSPEC 244 HCOMP(J,I) = MAX(HCOMP(J,I),0.0D0) 245 HCOMP(J,I) = MIN(HCOMP(J,I),1.0D0) 246 175 175 CONTINUE 247 HCOMP(3,I) = MIN(HCOMP(3,I),1.0D0-HCOMP(1,I)-HCOMP(4,I)) 248 HCOMP(2,I) = 1.0D0 - HCOMP(1,I) - HCOMP(3,I) - HCOMP(4,I) 249 180 180 CONTINUE 250 RETURN 251 END

1"""Python version of rotmix.f (YREC 7be4133): nuclear burning and 2gravitational settling for one step of the rotational-mixing loop.""" 3 4import math 5 6import numpy as np 7 8from pyrec import constants_mod 9from pyrec import dimensions_mod 10from pyrec import iounits_mod 11from pyrec import mdphy_mod 12from pyrec import nmlpar_mod 13from pyrec import scrtch_mod 14from pyrec.grsett import grsett 15from pyrec.kemcom import kemcom 16from pyrec.microdiff import microdiff 17from pyrec.rates import rates 18 19JSON = dimensions_mod.json 20 21# Local arrays (SAVE): the rates from rates, the settling's 22# pressure gradient and flags, and the saved SDEL(2,.). 23hr1 = np.zeros(JSON + 1) 24hr2 = np.zeros(JSON + 1) 25hr3 = np.zeros(JSON + 1) 26hr4 = np.zeros(JSON + 1) 27hr5 = np.zeros(JSON + 1) 28hr6 = np.zeros(JSON + 1) 29hr7 = np.zeros(JSON + 1) 30hr8 = np.zeros(JSON + 1) 31hr9 = np.zeros(JSON + 1) 32hr10 = np.zeros(JSON + 1) 33hr11 = np.zeros(JSON + 1) 34hr12 = np.zeros(JSON + 1) 35hr13 = np.zeros(JSON + 1) 36hf1 = np.zeros(JSON + 1) 37hf2 = np.zeros(JSON + 1) 38hqpr = np.zeros(JSON + 1) 39temp = np.zeros(JSON + 1) 40lcz = np.zeros(JSON + 1, dtype=bool) 41 42 43def rotmix(delts: float, hcomp: np.ndarray, hs2: np.ndarray, ht: np.ndarray, m: int, mrzone: np.ndarray, 44 mxzone: np.ndarray, nrzone: int, nzone: int, hstot: float, hd: np.ndarray, hs: np.ndarray, hr: np.ndarray, 45 hp: np.ndarray, lc: np.ndarray, hs1: np.ndarray) -> None: 46 """Burn each hot radiative shell and each convection zone with kemcom 47 (which also homogenizes the zones), apply the settling of helium 48 (microdiff or grsett, in substeps limited by the settling time at the 49 outer point), and clip and renormalize the composition. 50 51 Translates: rotmix.f#ROTMIX (7be4133) 52 53 Arguments: the step delts (s); the composition hcomp[species, shell] 54 (changed in place); shell masses hs2; log T ht; m shells; the zone 55 lists mrzone[k, 1..2] (radiative) and mxzone[k, 1..2] (convective), 56 with their counts nrzone, nzone; for the settling, log total mass 57 hstot, hd, hs, hr, hp (logs), the convection flags lc and the masses 58 hs1 (g). None is set (nmlpar_mod.ldify may be switched off). 59 60 """ 61 62 c = constants_mod 63 n = nmlpar_mod 64 65 # NSPEC is the number of species being tracked (Fortran lines 74–78). 66 if n.lexcom: 67 nspec = 15 68 else: 69 nspec = 11 70 71 # The step in years (Fortran line 80). 72 ddage = delts / c.csecyr 73 74 # The burning rates on the model passed in (hd, ht and the sub-step's 75 # own composition), as mix does (MHP via Claude, 10/29/26) (Fortran lines 87–110). 76 for ii in range(1, m + 1): 77 if ht[ii] <= n.tcut[1]: 78 for arr in (hr1, hr2, hr3, hr4, hr5, hr6, hr7, hr8, hr9, hr10, 79 hr11, hr12, hr13, hf1, hf2): 80 arr[ii] = 0.0 81 else: 82 rates(hd[ii], ht[ii], hcomp[1, ii], hcomp[2, ii], hcomp[4, ii], 83 hcomp[5, ii], hcomp[6, ii], hcomp[7, ii], hcomp[9, ii], 84 hcomp[11, ii], ii, hr1, hr2, hr3, hr4, hr5, hr6, hr7, 85 hr8, hr9, hr10, hr11, hr12, hr13, hf1, hf2) 86 87 # Radiative zones, shell by shell; stop at the first shell below the 88 # nuclear cutoff (GOTO 45 leaves both loops) (Fortran lines 118–130). 89 # Question (see the page): any later radiative zone is skipped. 90 cool = False 91 for k in range(1, nrzone + 1): 92 for j in range(mrzone[k, 1], mrzone[k, 2] + 1): 93 if ht[j] <= n.tcut[1]: 94 cool = True 95 break 96 ibegin = j 97 iend = j 98 kemcom(ht, ibegin, iend, hr1, hr2, hr3, hr4, hr5, hr6, hr7, 99 hr8, hr10, hr11, hr12, hf1, hs2, hcomp, ddage) 100 if cool: 101 break 102 103 # label 45. Convection zones, each as a unit (kemcom homogenizes it) (Fortran lines 135–142). 104 for k in range(1, nzone + 1): 105 ibegin = mxzone[k, 1] 106 iend = mxzone[k, 2] 107 kemcom(ht, ibegin, iend, hr1, hr2, hr3, hr4, hr5, hr6, hr7, 108 hr8, hr10, hr11, hr12, hf1, hs2, hcomp, ddage) 109 110 # Microscopic diffusion of helium (Fortran lines 151–237). 111 if n.ldify: 112 _settle(delts, hcomp, ht, m, mrzone, mxzone, nrzone, nzone, 113 hstot, hd, hs, hr, hp, lc, hs1) 114 # End of gravitational settling (label 170). 115 116 # Renormalize the composition (small negative abundances...) (Fortran lines 242–249). 117 for i in range(1, m + 1): 118 for j in range(1, nspec + 1): 119 hcomp[j, i] = max(hcomp[j, i], 0.0) 120 hcomp[j, i] = min(hcomp[j, i], 1.0) 121 hcomp[3, i] = min(hcomp[3, i], 1.0 - hcomp[1, i] - hcomp[4, i]) 122 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 123 124 125def _settle(delts: float, hcomp: np.ndarray, ht: np.ndarray, m: int, mrzone: np.ndarray, mxzone: np.ndarray, 126 nrzone: int, nzone: int, hstot: float, hd: np.ndarray, hs: np.ndarray, hr: np.ndarray, hp: np.ndarray, 127 lc: np.ndarray, hs1: np.ndarray) -> None: 128 """The settling block inside IF(LDIFY); each early return is a 129 GOTO 170. 130 131 """ 132 133 c = constants_mod 134 n = nmlpar_mod 135 sdel = scrtch_mod.sdel 136 137 # Settling stops for good once the central X is below XMIN (Fortran lines 152–155). 138 if hcomp[1, 1] < n.xmin: 139 n.ldify = False 140 return 141 142 # Locate the outer boundary (Fortran lines 163–180). 143 if not lc[m]: 144 text = ("NO SURFACE CZ - DIFFUSION NOT MEANINGFUL\n" 145 "STOPPED IN SUBROUTINE MIX") # format 911 146 # Suspected bug (see the page): the message names MIX. 147 print(" " + text) 148 iounits_mod.write(iounits_mod.ishort, " " + text) 149 raise SystemExit(1) # STOP 1, the run-time error code 150 # Exit if overshoot makes the star fully convective (MHP 7/92). 151 if mxzone[nzone, 1] <= 2 and mxzone[nzone, 2] == m: 152 return 153 jmax = mrzone[nrzone, 2] 154 for j in range(jmax, 0, -1): 155 if hcomp[2, j] > n.ymin: 156 break # go to 150 157 else: 158 return # Y < YMIN for the whole star 159 # label 150. Suspected bug (see the page): the loop above moves J, not 160 # JMAX, so JMAX stays the top of the last radiative zone and the first 161 # shell with Y > YMIN is never used. 162 163 # The pressure gradient; SDEL(2,.) holds DELM meanwhile (Fortran lines 181–186). 164 stot = math.exp(c.cln * hstot) 165 for i in range(1, m + 1): 166 temp[i] = sdel[2, i] 167 sdel[2, i] = mdphy_mod.delm[i] 168 hqpr[i] = -math.exp(c.cln * (hd[i] + c.cgl + hs[i] 169 - 2.0 * hr[i] - hp[i])) 170 171 # Only check for timestep cutting if JMAX is large (BC 6/92) (Fortran lines 189–213). 172 if jmax >= 2: 173 # FM is the mass fraction above the outer point. 174 fm = (stot - hs1[jmax]) / stot 175 # Settling time of helium there (Michaud et al. 1984). 176 tscale = 4.348e21 * c.csecyr * fm / math.exp(c.cln * 1.5 * ht[jmax]) 177 dtmax = n.dt_gs * tscale 178 nstep = int(delts / dtmax) 179 if math.fmod(delts, dtmax) != 0.0 or nstep == 0: 180 nstep = nstep + 1 181 dt = delts / float(nstep) 182 else: 183 nstep = 1 184 dt = delts 185 186 # The convection flags from the zone lists (Fortran lines 214–223). 187 for j in range(1, nrzone + 1): 188 for i in range(mrzone[j, 1], mrzone[j, 2] + 1): 189 lcz[i] = False 190 for j in range(1, nzone + 1): 191 for i in range(mxzone[j, 1], mxzone[j, 2] + 1): 192 lcz[i] = True 193 194 # The settling substeps (Fortran lines 224–232). 195 for ns in range(1, nstep + 1): 196 # LNEWDIF: the new routines in microdiff; else grsett. 197 if n.lnewdif: 198 # Since 7be4133 microdiff and grsett work on copies of dt 199 # and stot and leave them unchanged. 200 microdiff(dt, hcomp, hqpr, hr, hd, hs1, ht, lcz, m, stot) 201 else: 202 grsett(dt, hcomp, hqpr, hr, hd, hs1, ht, lcz, m, stot) 203 204 # SDEL(2,.) restored (Fortran lines 233–235). 205 for i in range(1, m + 1): 206 sdel[2, i] = temp[i]