crrect.f

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

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

Variables on this page (22)

cln, fcorr, fcorr0, fcorri, hcompp, htoler, iowr, ishort, json, lcorr, lczo, ldh, lexcom, lkuthe, lrot, senv, xnew, xxdh, yydh, zdh, znew, zzdh

Purpose

This file runs the Henyey iterations that correct the interior model at one level of the time step: it sets the surface boundary condition from the envelope, mixes and burns the composition, builds the linearized structure equations, solves them for the corrections to $\log P$, $\log T$, $\log r$ and $L$, applies the corrections (damped), and decides whether the model has converged.

How it fits

Inputs: the arguments: the model (HP, HT, HR, HL, HD, HS, HS1, HS2, HCOMP, HPP, HTT), DELTS (time step, s), NITER (the maximum number of iterations), ITLVL (the level of iteration within the step), the flags LNEW, LRESET, LSBC, the envelope triangle (envelopes integrated at three points around the model's $(\log T_\mathrm{eff}, \log L)$, whose values at the fitting point are fitted by a plane to give the outer boundary condition; SURFBC: TRIT, TRIL, PS, TS, RS, FTRI, ISTORE, ESTORE), rotation (FP, FT, OMEGA, HJM, …), and CMIXL (the mixing length, set by main.f from the namelist's CMIXLA, passed on to SURFBC, MIX and COEFFT); from NMLPAR_MOD, the tolerances HTOLER, the damping settings FCORR0, FCORRI, and the options LKUTHE (the helium flash), LCORR (print corrections), LROT, LDH, LEXCOM; the current damping FCORR (CTOL_MOD); the start-of-step composition HCOMPP and, with rotation, the convective flags LCZO that HPOINT (or STARIN, for the starting model) stored from OVROT (OLDMOD_MOD), XNEW, ZNEW (COMP_MOD), the Debye–Hückel composition (DEBHU_MOD), and the output units from IOUNITS_MOD.
Who calls it: the main program (main.f), within the evolution loop.
What it calls: SURFBC (the surface boundary condition from envelope integrations, as fitting coefficients CFENV), MIX (burning and mixing over the step), COEFFT (the equations and their derivatives, reduced by REDUCE), HSOLVE (the back-substitution), and with rotation GETROT and FPFT.
Outputs: the corrected model, BL, TEFFL, LFINI (converged), LARGE (corrections too large: main.f then calls STARIN, which stops the run; see main.f, step 2 of the evolution step), ITDONE, HMAX/NHMAX (largest corrections and where), TLUMX (luminosity by source, renormalized), and with rotation FP, FT, OMEGA, HKEROT. With LCORR, a line per iteration to ISHORT.

The method

The system as a whole. The unknowns are the corrections $\delta\log P$, $\delta\log T$, $\delta\log r$ and $\delta L$ at each of the $M$ points, $4M$ in all; the composition is not among them (MIX sets it: main.f's call before levels 1 and 2, and step 2 at levels 3 and 4). They are fixed by $4(M-1)$ difference equations, one for each structure equation between each pair of neighbouring points; two central conditions at point 1 (COEFFT); and two surface conditions, the envelope fit below (HSRF). The matrix is block-bidiagonal, and REDUCE and HSOLVE are the forward and backward sweeps that solve it.

If LSBC, call SURFBC with the envelope composition. It returns the fit of the envelope's base as linear relations in $\log P$ and $\log T$ at the outermost interior point: $$\log r = c_1\log P + c_2\log T + c_3, \qquad \log L = c_4\log P + c_5\log T + c_6, \qquad \log T_\mathrm{eff} = c_7\log P + c_8\log T + c_9.$$ With LDH, the composition for the Debye–Hückel correction (XXDH, YYDH, ZZDH, ZDH) is first set from shell M.
At ITLVL > 2 with $\Delta t > 0$: reset the composition to its start-of-step value and call MIX. The call to OVROT that followed is commented out (lines 162–165 ↗); with rotation, OVROT runs in HPOINT after each rezone ("CALCULATE OVERSHOOT," hpoint.f, lines 779–781 ↗) and in STARIN for the starting model (starin.f, lines 1032–1033 ↗), and CRRECT uses the flags they store, LCZO (settled from the code, Nov 9, 2026).
If the boundary condition was recomputed, the triangle did not flip (LFLP from SURFBC) and the model had converged at the previous level (LFINI), return at once.
Iterate (at most NITER times):
COEFFT and REDUCE build the system. Renormalize TLUMX so that its sum equals $L$ at the surface (not during the helium flash).
From the second iteration, if the largest residuals of the equations (HMAX from REDUCE) are below HTOLER(5,·), the model has converged.

Question 1 (lines 190–201 ↗). Residual test before solving. At iteration 2 and later, if the residuals are below HTOLER(5,·), the routine returns as converged without solving. The comment notes that boundary equations are not checked. Can the surface condition then be unmet? (The message prints ITDONE+1 as the iteration count, but ITDONE is not incremented on this path.)
The surface rows from the fit above: HSRF(1..3) linearize the $\log r$ relation and HSRF(4..6) the $\log L$ relation, multiplied by $\ln 10\,L$ since $L$, not $\log L$, is the unknown; the $T_\mathrm{eff}$ relation only updates TEFFL afterwards. HSOLVE gives the corrections $\delta$.
Converged if $\max|\delta\log P|$, $\max|\delta\log T|$, $\max|\delta\log r|$ and $\max \min(|\delta L|, |\delta L/L|)$ are all below HTOLER(·,1); stop (return with LARGE) if any exceeds HTOLER(·,2).

Question 2 (line 234 ↗). The luminosity test. The convergence test on $L$ uses $\min(|\delta L|, |\delta L/L|)$, with HL in solar units: where $L > 1\,L_\odot$ this is the relative change, elsewhere the absolute one (in a star fainter than the Sun, the absolute one everywhere); the first shell enters with $|\delta L|$ alone (line 215 ↗). Is the mixed test intended?
Apply $x \leftarrow x + f\,\delta x$, with the damping $f$ = FCORR, raised by FCORRI per iteration up to 1 if FCORR0 > 0; during the helium flash $f = 0.8$ once $|\delta L| \le 0.5$ and 1 once $\le 5\times10^{-3}$. The changes of $\log P$, $\log T$ over the step (HPP, HTT) are updated too.
New $\log L$ and $\log T_\mathrm{eff}$ (from the fit); with rotation, new $\Omega$, $f_P$, $f_T$ and rotational kinetic energy $\tfrac12\Omega j\,\Delta m$. GETROT gets the flags LCZO, so that it imposes the convection-zone rotation law on each convection zone as a whole (WCZIMP) and inverts the other shells one by one.

From YREC's documentation

Physics Options (Numerical convergence options):

variable	default	meaning (documentation)
`HTOLER(1..4,1)`	6.0D-5, 4.5D-5, 3.0D-5, 9.0D-5	convergence tolerances for $\log P$, $\log T$, $\log R$, $L/L_\mathrm{tot}$: the model has converged when the largest correction anywhere is below all four (all local)
`HTOLER(5,1)`, `HTOLER(5,2)`	3.0D-5, 2.5D-6	if analytic and numerical derivatives agree better than this ($P$ and $R$; $T$ and $L$), the iteration is skipped
`HTOLER(1..4,2)`	0.9, 0.5, 0.5, 2.0	abort thresholds: a correction larger than this is a divergence (rarely used)
`FCORR0`	0.8	the predicted corrections can be overestimated at first; the first iteration multiplies them by this
`FCORRI`	0.1	the increase of that factor at each iteration (with the defaults, 0.8, 0.9, 1.0); larger is faster but less stable
`LCORR`	true	write the details of convergence to the `.short` file

The code does not compare derivatives for HTOLER(5,·): it compares the largest residuals of the difference equations (HMAX, from REDUCE) with HTOLER(5,1) for $P$ and $r$ and HTOLER(5,2) for $T$ and $L$, and returns as converged ("R.H.S. BELOW TOLERANCES") without solving (step 4).

HTOLER(5,·) bounds residuals, not derivatives. The Physics namelist page says that "If the difference between analytic and numerical derivatives is better than this, iteration is skipped." The code ("CHECK ON SIGNIFICANCE OF R.H.S. EQUATIONS," message "R.H.S. BELOW TOLERANCES": lines 190–196 ↗) compares with them HMAX, the largest residuals of the difference equations, which REDUCE fills (reduce.f, lines 33–40 ↗), HTOLER(5,1) for $P$ and $r$ and HTOLER(5,2) for $T$ and $L$, before HSOLVE's call (line 212 ↗) reuses HMAX for the corrections; no numerical derivatives are formed, so the documentation could describe the residual test.

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) `JENV` is unset at levels 1 and 2.

(line 57 ↗; lines 150–161 ↗; lines 176–179 ↗) JENV, the first shell of the convective envelope, is a local of CRRECT, set only by its own call to MIX at levels 3 and 4 (ITLVL > 2), and passed to COEFFT. COEFFT reads it when the entropy terms are on and mass is gained (LMDOT and DMDT0 > 0): shells with IM.GE.JENV get the extra changes DLOGT, DLOGP in $\log T$, $\log P$ (coefft.f, lines 324–331 ↗). At levels 1 and 2 of the first step, JENV has never been set; at levels 1 and 2 of later steps it is the value from the previous step's MIX, not the one the main program's own MIX and FINDSH have just computed (main.f keeps its own JENV). What follows: with accretion, the extra terms go to the wrong shells (or, with an unset value, to all or none) until level 3. Only with LMDOT (off by default) and DMDT0 > 0 (the default is $-10^{-14}$, mass loss).

A possible fix: pass the main program's JENV (set by STARIN, MIX and FINDSH) to CRRECT as an argument, in place of the local, so that CRRECT's MIX at levels 3 and 4 updates the same variable.

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 a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning, and 8 made-up ones; one of these, with FCORE = 1.0, stopped in its first second). Two made-up runs had mass loss on (LMDOT and LREIMER: YREC's 1 $M_\odot$ case to helium ignition, and a 3 $M_\odot$ main-sequence case), but mass loss makes the rate negative, and the path needs LMDOT with a mass gain (DMDT0 > 0, LREIMER off), which no case or document sets; the namelists' comment on LMDOT reads "Constant Mdot (negative) and Reimers are currently enabled." A run with accretion would reach it; none was run, and whether level 3 removes the effect is still not traced. No fix was run.

Definitions

subroutineCRRECT

lines 11–313 · GitHub

The Henyey iteration for one level of the time step: surface boundary condition (SURFBC), mixing and burning (MIX), equations (COEFFT), solution (HSOLVE), damped corrections and convergence tests against HTOLER.

Uses modules COMP_MOD, CONSTANTS_MOD, CTOL_MOD, DEBHU_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD, OLDMOD_MOD

Called by (file level) ×4

Flowchart of its calls, in order (6)

if (LSBC), lines 132–148

SURFBC 144

if (ITLVL.GT.2 .AND. DELTS.GT.0.0D0), lines 150–161

MIX 158

for ITER = 1,NITER (DO 100), lines 170–311

COEFFT 176→HSOLVE 212

if (LROT), lines 294–307

GETROT 301→FPFT 302

Change log

The changes since rp_1

Changed since rp_1. The new argument CMIXL (the mixing length, formerly in /CONST3/), passed on to SURFBC, MIX and COEFFT; COMMON blocks to modules, IMPLICIT NONE.

Changed since a8a5950. GETROT gets LCZO, the convective flags from OVROT, instead of the argument LCZ, which main.f never set, so that every shell was inverted on its own and the convection-zone law (WCZIMP) was never imposed in the structure solve; the argument LCZ, then unused, is removed, and so are QDT and QDP, now locals of COEFFT; the test that sets LSBC is parenthesized as the compiler parses it, LNEW .OR. (LRESET .AND. ITLVL.EQ.2), since a new run needs a triangle at level 1 (this settles former question 1). The later questions are renumbered.

1C JANUARY 28, 1991: 2C REPLACED EGEN SO THAT IT FOLLOWS FULL CNO NOT CN IN EQUILIBRIUM AND 3C FEBRUARY 11, 1991: 4C REPLACED CONVEC SO THAT WOULD CALCULATE OVERSHOOT CORRECTLY. 5C ADDED ROUTINES OVERSH AND HSUBP. 6C 7C 8C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 9C CRRECT 10C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 11 SUBROUTINE CRRECT(DELTS,M,NITER,LFINI,LARGE,LNEW,LRESET,LSBC, 12 * *TRIT,TRIL,CFENV,PS,TS,RS,FTRI,ISTORE,ESTORE,HSTOT,BL,TEFFL, 13 * *HD,HHA,HHB,HHC,HL,HMAX,HP,HPP,HR,HS,HS1,HS2,HSRF,HT,HTT,HCOMP,LC, 14C *NHMAX,TLUMX,LATMO,LDERIV,LMIX,LOCOND,QDT,QDP,MODEL,FP,FT,ETA2, 15C *OMEGA,R0,ITDONE,HG,HI,HJM, ITLVL,LCZ,MRZONE,MXZONE,NRZONE,NZONE, ! KC 2025-05-31 16 * *NHMAX,TLUMX,LATMO,LDERIV,LMIX,LOCOND,FP,FT,ETA2, 17 * *OMEGA,R0,ITDONE,HG,HI,HJM, ITLVL,MXZONE,QIW,HKEROT,HKEROT0, 18 * *CMIXL) 19C 10/30/26 MHP via Claude 20C code changed to drop the unused LCZ argument. 21C 11/02/26 MHP via Claude 22C code changed to drop the QDT and QDP arguments, now COEFFT locals. 23C 10/15/26 MHP via Claude 24C code changed to take JSON from DIMENSIONS_MOD instead of declaring 25C PARAMETER(JSON=5000) here. 26 USE DIMENSIONS_MOD 27C 10/14/26 MHP via Claude 28C code changed to convert COMMON/OLDMOD/ into OLDMOD_MOD. 29 USE OLDMOD_MOD 30C 10/14/26 MHP via Claude 31C code changed to convert COMMON/DEBHU/ into DEBHU_MOD. 32 USE DEBHU_MOD 33C 10/14/26 MHP via Claude 34C code changed to convert COMMON/CTOL/ into CTOL_MOD. 35 USE CTOL_MOD 36C 10/14/26 MHP via Claude 37C code changed to convert COMMON/COMP/ into COMP_MOD. 38 USE COMP_MOD 39C 10/14/26 MHP via Claude 40C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 41 USE NMLPAR_MOD 42C 10/14/26 MHP via Claude 43C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 44 USE CONSTANTS_MOD 45C 10/14/26 MHP via Claude 46C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 47 USE IOUNITS_MOD 48C 10/20/26 MHP via Claude 49C code changed to replace implicit typing with IMPLICIT NONE and 50C explicit declarations. 51 IMPLICIT NONE 52 REAL*8 BL,CFENV,CMIXL,DDAGE,DELTS,ESTORE,ETA2,FACT,FP,FPL,FT,FTL, 53 * * FTRI,HCOMP,HD,HG,HHA,HHB,HHC,HI,HJ,HJM,HKEROT,HKEROT0,HL, 54 * * HLMAX,HMAX,HP,HPP,HR,HS,HS1,HS2,HSRF,HSTOT,HT,HTT,OMEGA, 55 * * PLIM,PS,QIW,R0,RS,TEFFL,TEMP,TEST,THL,TLUMH,TLUMHE, 56 * * TLUMX,TRIL,TRIT,TS,X,Z 57 INTEGER I,IJ,ISTORE,ITDONE,ITER,ITLVL,J,JCORE,JENV,KATM,KENV, 58 * * KSAHA,M,MXZON0,MXZONE,NHMAX,NITER,NSPEC 59 LOGICAL*4 LARGE,LATMO,LC,LDERIV,LFINI,LFLP,LMIX,LNEW,LOCOND, 60 * * LRESET,LSBC 61 62C INPUTS ASSUMES GIVEN LOG(TE) AS TEFFL 63C ASSUMES GIVEN LOG(L/LSUN) AS BL 64C DELTS = TIME STEP IN SECONDS 65C NITER = ITERATION NUMBER 66C OUTPUTS LFINI = .T. IF MODEL HAS CONVERGED 67C LARGE = .T. IF CORRECTIONS ARE TOO LARGE 68C 9/25/26 MHP via Claude 69C code changed to consolidate all I/O logical unit numbers, previously 70C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 71C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 72C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 73C COMMON/IOUNITS/ 74C 10/11/26 MHP via Claude 75C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared 76C positionally in all 17 files sharing this block but never read or 77C written in executable code anywhere, and no namelist file sets it. 78C 10/12/26 MHP via Claude 79C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 80C namelist variables that the physics routines read from 27 COMMON 81C blocks which also held derived or working members. Those blocks keep 82C only their non-namelist members; five of them -- NEWMX, VMULT, 83C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 84C Derivation-only namelist variables from the same blocks went to 85C COMMON/NMLBYPASS/ instead, not here. 86C 87C Declared identically in all 121 files that hold any of these members. 88C Generated from one member list, so every copy is byte-identical; the 89C linker's 'size of symbol changed' warning is the check that they are. 90C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration 91C (no member referenced here; residue of the obsolete-option removals). 92C 10/13/26 MHP via Claude: removed the unused COMMON/CENV/ declaration 93C (no member referenced here; residue of the obsolete-option removals). 94C 9/25/26 MHP via Claude 95C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 96C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 97C DBG 7/92 COMMON BLOCK ADDED TO COMPUTE DEBYE-HUCKEL CORRECTION. 98 DIMENSION MXZONE(12,2),MXZON0(12,2) 99 DIMENSION TRIT(3),TRIL(3),CFENV(9),PS(3),TS(3),RS(3),ESTORE(4), 100 * *TLUMX(8),FP(JSON),FT(JSON),ETA2(JSON),OMEGA(JSON),R0(JSON) 101 DIMENSION HD(JSON),HHA(4,2,JSON),HHB(4,JSON),HHC(JSON),HL(JSON), 102 * * HMAX(4),HP(JSON),HPP(JSON),HR(JSON),HS(JSON),HS1(JSON),HS2(JSON), 103 * * HSRF(6),HT(JSON),HTT(JSON),HCOMP(15,JSON),LC(JSON),NHMAX(4), 104 * * HG(JSON),HI(JSON),HJM(JSON) 105C MHP 10/02 ADDED MISSING DIMENSION STATEMENTS 106 DIMENSION HKEROT(JSON),HKEROT0(JSON), 107C * LCZ(JSON),QIW(JSON),MRZONE(13,2) ! KC 2025-05-31 108 * * QIW(JSON) 109 SAVE 110 111 IF(NITER.LE.0) RETURN 112 BL = DLOG10(HL(M)) 113C ZERO COUNTERS 114 KENV = 0 115 KATM = 0 116 KSAHA = 0 117 SENV = HS(M) - HSTOT 118C 10/30/26 MHP via Claude 119C code changed to parenthesize as the compiler parses it. A new run 120C needs a triangle at level 1; MAIN already sets LSBC at level 2. 121 IF(LNEW .OR. (LRESET .AND. ITLVL.EQ.2)) LSBC = .TRUE. 122C FIND NEW FP AND FT IF MODEL IS ROTATING 123 IF(LROT.AND.LSBC) THEN 124 FPL = FP(M) 125 FTL = FT(M) 126 ELSE 127 FPL = 1.0D0 128 FTL = 1.0D0 129 ENDIF 130C SET UP SURFACE BOUNDARY CONDITIONS-2ND AND 3RD LEVELS OF ITER ONLY 131C FIND ENVELOPE MASS AND SET X AND Z TO ENVELOPE VALUES 132 IF(LSBC) THEN 133 X = XNEW 134 Z = ZNEW 135 PLIM = HP(M) 136 IF (LDH) THEN 137 XXDH = HCOMP(1,M) 138 YYDH = HCOMP(2,M)+HCOMP(4,M) 139 ZZDH = HCOMP(3,M) 140 ZDH(1) = HCOMP(5,M)+HCOMP(6,M) 141 ZDH(2) = HCOMP(7,M)+HCOMP(8,M) 142 ZDH(3) = HCOMP(9,M)+HCOMP(10,M)+HCOMP(11,M) 143 END IF 144 CALL SURFBC(TRIT,TRIL,CFENV,PS,TS,RS,FTRI,ISTORE,ESTORE, 145C * LNEW,LRESET,LSBC,KSAHA,KENV,KATM,HSTOT,BL, ! KC 2025-05-31 146 * * LNEW,LRESET,KSAHA,KENV,KATM,HSTOT,BL, ! LSBC 147 * * TEFFL,X,Z,FPL,FTL,LFLP,PLIM,LC,M,CMIXL) 148 ENDIF 149C 7/91 ADD CALL TO MIX 150 IF(ITLVL.GT.2 .AND. DELTS.GT.0.0D0)THEN 151 NSPEC=11 152 IF(LEXCOM)NSPEC=15 153 DO 2 I = 1,M 154 DO 1 J = 1,NSPEC 155 HCOMP(J,I) = HCOMPP(J,I) 156 1 1 CONTINUE 157 2 2 CONTINUE 158 CALL MIX(DELTS,HCOMP,HD,HL,HP,HR,HS,HS1,HS2,HSTOT, 159C * HT,ITLVL,LC,M,QDT,QDP,DDAGE,JCORE,JENV, ! KC 2025-05-31 160 * * HT,ITLVL,LC,M,DDAGE,JCORE,JENV,MXZONE,MXZON0,TEFFL,CMIXL) 161 ENDIF 162c IF(LROT)THEN 163c CALL OVROT(HCOMP,HD,HP,HR,HS,HT,LC,M,LCZ,MRZONE,MXZONE, 164c * NRZONE,NZONE) 165c ENDIF 166C IF MODEL CONVERGED ON THE PREVIOUS LEVEL OF ITERATION AND 167C THE TRIANGLE WAS CHECKED AND NOT FLIPPED,SKIP CORRECTION ROUTINE. 168 IF(LSBC .AND. .NOT.LFLP .AND. LFINI) RETURN 169 LFINI = .FALSE. 170 DO 100 ITER = 1,NITER 171C DO HENYEY REDUCTION 172 HMAX(1) = 0.0D0 173 HMAX(2) = 0.0D0 174 HMAX(3) = 0.0D0 175 HMAX(4) = 0.0D0 176 CALL COEFFT(DELTS,M,HD,HHA,HHB,HHC,HL,HMAX,HP,HPP,HR,HS,HS1, 177 * * HS2,HT,HTT,HCOMP,LC,TLUMX,LATMO,LDERIV,LMIX,LOCOND, 178C * KSAHA,MODEL,FP,FT,HKEROT,HKEROT0,JENV,TEFFL) ! KC 2025-05-31 179 * * KSAHA,FP,FT,HKEROT,HKEROT0,JENV,TEFFL,CMIXL) 180C RENORMALIZE TLUMX-S 181CCC TAKE OUT RENORMALIZATION DURING HE FLASH (NON-THERMAL EQUALIBRIUM) 182 THL = TLUMX(1)+TLUMX(2)+TLUMX(3)+TLUMX(4)+TLUMX(5)+TLUMX(6) 183 * * +TLUMX(7)+TLUMX(8) 184 IF (.NOT.LKUTHE .AND. THL.GT.0.0D0) THEN 185 TEMP = HL(M)/THL 186 DO 10 J = 1,8 187 TLUMX(J) = TLUMX(J)*TEMP 188 10 10 CONTINUE 189 ENDIF 190C CHECK ON SIGNIFICANCE OF R.H.S. EQUATIONS FOR P AND T 191C N.B. DOES NOT CHECK DIFFERENCES IN BOUNDARY EQUATIONS 192 IF(ITER.GT.1) THEN 193 IF(HMAX(1).LE.HTOLER(5,1) .AND. HMAX(2).LE.HTOLER(5,2) .AND. 194 * * HMAX(3).LE.HTOLER(5,1) .AND. HMAX(4).LE. HTOLER(5,2)) THEN 195 WRITE(ISHORT,20) (HMAX(J),J=1,4) 196 20 20 FORMAT(' R.H.S. BELOW TOLERANCES--P',1PE9.2,' T ',E9.2, 197 * * ' R ',E9.2,' L ',E9.2) 198 WRITE(ISHORT,75) ITDONE+1 199 LFINI = .TRUE. 200 RETURN 201 ENDIF 202 ENDIF 203C SET UP HENYEY COEFFICIENTS FOR SURFACE 204 HSRF(1) = -CFENV(1) 205 HSRF(2) = -CFENV(2) 206 HSRF(3) = CFENV(1)*HP(M) + CFENV(2)*HT(M) + CFENV(3) - HR(M) 207 TEMP = CLN*HL(M) 208 HSRF(4) =-TEMP*CFENV(4) 209 HSRF(5) =-TEMP*CFENV(5) 210 HSRF(6) = TEMP*(CFENV(4)*HP(M)+CFENV(5)*HT(M)+CFENV(6)-BL) 211C DO BACK SOLUTION FOR CORRECTIONS 212 CALL HSOLVE(M,HHA,HHB,HSRF) 213C CHECK ON MAXIMUM CORRECTIONS 214 DO 30 J = 1,4 215 HMAX(J) = DABS(HHB(J,1)) 216 NHMAX(J) = 1 217 30 30 CONTINUE 218 DO 40 I = 2,M 219 TEST = DABS(HHB(1,I)) 220 IF(HMAX(1).LE.TEST) THEN 221 HMAX(1) = TEST 222 NHMAX(1) = I 223 ENDIF 224 TEST = DABS(HHB(2,I)) 225 IF(HMAX(2).LE.TEST) THEN 226 HMAX(2) = TEST 227 NHMAX(2) = I 228 ENDIF 229 TEST = DABS(HHB(3,I)) 230 IF(HMAX(3).LE.TEST) THEN 231 HMAX(3) = TEST 232 NHMAX(3) = I 233 ENDIF 234 TEST = DMIN1(DABS(HHB(4,I)),DABS(HHB(4,I)/HL(I))) 235 IF(HMAX(4).LE.TEST) THEN 236 HMAX(4) = TEST 237 NHMAX(4) = I 238 ENDIF 239 40 40 CONTINUE 240CCC HE FLASH -- OK FOR ALL 241 HLMAX = HMAX(4) 242C LFINI = T IF MAX CORRECTIONS LESS THAN CONVERGENCE CRITERIA SET IN 243C HTOLER. LARGE = T IF MAX CORRECTIONS GREATER THAN LARGEST CORRECTIONS 244C ALLOWED, ALSO SET IN HTOLER 245 LFINI = HMAX(1).LT.HTOLER(1,1) .AND. HMAX(2).LT.HTOLER(2,1) 246 * * .AND. HMAX(3).LT.HTOLER(3,1) .AND. HMAX(4).LT.HTOLER(4,1) 247 LARGE = HMAX(1).GT.HTOLER(1,2) .OR. HMAX(2).GT.HTOLER(2,2) 248 * * .OR. HMAX(3).GT.HTOLER(3,2) .OR. HMAX(4).GT.HTOLER(4,2) 249 DO 50 J = 1,4 250 IJ = NHMAX(J) 251 HMAX(J) = HHB(J,IJ) 252 50 50 CONTINUE 253 IF(FCORR0.GT.0.0D0) FCORR = DMIN1(1.D0,FCORR+FCORRI) 254C HE FLASH CHANGE 255 FACT = FCORR 256 IF (LKUTHE) THEN 257 IF (HLMAX.LE.5.0D-1) FACT=8.0D-1 258 IF (HLMAX.LE.5.0D-3) FACT=1.0D0 259 ENDIF 260 TLUMH = TLUMX(1) + TLUMX(2) + TLUMX(3) + TLUMX(4) 261 TLUMHE= TLUMX(5) + TLUMX(8) 262 IF(LCORR) THEN 263 WRITE (ISHORT,60) LFINI,HMAX(4),HTOLER(4,1),NHMAX(4) 264 60 60 FORMAT (1X,'DEL-L/L ',L2,1P2E12.4,5X,I5) 265 WRITE(ISHORT,70)(NHMAX(J),HMAX(J),J=1,4),FACT,TLUMH, 266 * * TLUMHE,(TLUMX(J),J=6,7) 267 70 70 FORMAT(' CORR',I5,'P',1PE9.2,I5,'T',E9.2,I5,'R',E9.2,I5,'L', 268 * * E9.2,' F=',0PF5.3,' E-HY',1PE10.3,' HE',E10.3,' NU',E10.3, 269 * * ' G',E10.3) 270 IF(LFINI) THEN 271 WRITE(ISHORT,75) ITDONE+1 272 75 75 FORMAT(10X,'MODEL CONVERGED AFTER ',I4,' ITERATIONS') 273 ENDIF 274 ENDIF 275 IF(LARGE) THEN 276 WRITE(IOWR,80) 277 WRITE(ISHORT,80) 278 80 80 FORMAT(1X,'-----CORRECTIONS EXCEEDED TOLERANCES') 279 RETURN 280 ENDIF 281C APPLY CORRECTIONS 282 DO 90 I = 1,M 283 TEMP = FACT*HHB(1,I) 284 HP(I) = HP(I) + TEMP 285 HPP(I) = HPP(I) + TEMP 286 TEMP = FACT*HHB(2,I) 287 HT(I) = HT(I) + TEMP 288 HTT(I) = HTT(I) + TEMP 289 HR(I) = HR(I) + FACT*HHB(3,I) 290 HL(I) = HL(I) + FACT*HHB(4,I) 291 90 90 CONTINUE 292 BL = DLOG10(HL(M)) 293 TEFFL = CFENV(7)*HP(M) + CFENV(8)*HT(M) + CFENV(9) 294 IF(LROT)THEN 295C 10/29/26 MHP via Claude 296C GETROT now gets LCZO, the rotation convective flags HPOINT (or STARIN, 297C for the starting model) stored from OVROT. It used to get LCZ, which 298C MAIN declares and passes but never sets, so every shell was inverted 299C on its own and WCZ's convection-zone law was never imposed here. 300C He-flash steps skip HPOINT and keep the previous model's flags. 301 CALL GETROT(HD,HJM,HR,HS,HS2,LCZO,M,ETA2,HI,OMEGA,QIW,R0) 302 CALL FPFT(HD,HR,HS,M,OMEGA,ETA2,FP,FT,HG,R0) 303 DO I = 1,M 304 HJ = HJM(I)*HS2(I) 305 HKEROT(I) = 0.5D0*OMEGA(I)*HJ 306 END DO 307 ENDIF 308 ITDONE = ITDONE + 1 309 IF(LFINI) RETURN 310 100 100 CONTINUE 311 312 RETURN 313 END

1"""Python version of crrect.f (YREC 7be4133): the Henyey iterations that 2correct the interior model at one level of the time step.""" 3 4import math 5 6import numpy as np 7 8from pyrec import comp_mod 9from pyrec import constants_mod 10from pyrec import ctol_mod 11from pyrec import debhu_mod 12from pyrec import iounits_mod 13from pyrec import nmlpar_mod 14from pyrec import oldmod_mod 15from pyrec.coefft import coefft 16from pyrec.fpft import fpft 17from pyrec.getrot import getrot 18from pyrec.hsolve import hsolve 19from pyrec.mix import mix 20from pyrec.surfbc import surfbc 21 22# CRRECT has SAVE. Locals kept between calls that matter: JCORE, JENV and 23# MXZON0 are set only by MIX (at ITLVL > 2), but JENV is passed to COEFFT at 24# every level; at levels 1 and 2 it keeps its value from an earlier call 25# (unset at the first call; 0 here). Suspected bug (see the page). 26# Fortran lines 52–60, 98–109. 27jcore = 0 28jenv = 0 29mxzon0 = np.zeros((13, 3), dtype=int) # MXZON0(12,2), index 0 unused 30ddage = 0.0 31lflp = False 32 33 34def _l(flag: bool) -> str: 35 """Fortran's L format of a logical.""" 36 37 return "T" if flag else "F" 38 39 40def crrect(delts: float, m: int, niter: int, lfini: bool, large: bool, lnew: bool, lreset: bool, lsbc: bool, 41 trit: np.ndarray, tril: np.ndarray, cfenv: np.ndarray, ps: np.ndarray, ts: np.ndarray, rs: np.ndarray, 42 ftri: float, istore: int, estore: np.ndarray, hstot: float, bl: float, teffl: float, hd: np.ndarray, 43 hha: np.ndarray, hhb: np.ndarray, hhc: np.ndarray, hl: np.ndarray, hmax: np.ndarray, hp: np.ndarray, 44 hpp: np.ndarray, hr: np.ndarray, hs: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, hsrf: np.ndarray, 45 ht: np.ndarray, htt: np.ndarray, hcomp: np.ndarray, lc: np.ndarray, nhmax: np.ndarray, tlumx: np.ndarray, 46 latmo: bool, lderiv: bool, lmix: bool, locond: bool, fp: np.ndarray, ft: np.ndarray, eta2: np.ndarray, 47 omega: np.ndarray, r0: np.ndarray, itdone: int, hg: np.ndarray, hi: np.ndarray, hjm: np.ndarray, itlvl: int, 48 mxzone: np.ndarray, qiw: np.ndarray, hkerot: np.ndarray, hkerot0: np.ndarray, cmixl: float) -> tuple[bool, 49 bool, bool, bool, bool, float, int, float, float, bool, bool, bool, bool, int]: 50 """Henyey iteration for one level of the time step: the surface 51 boundary condition (SURFBC), mixing and burning (MIX), the equations 52 (COEFFT), their solution (HSOLVE), damped corrections and the 53 convergence tests against HTOLER. 54 55 Translates: crrect.f#CRRECT (7be4133) 56 57 Arguments: delts, the time step (s); m, the number of shells; niter, 58 the maximum number of iterations; lfini (read: converged at the 59 previous level); the flags lnew, lreset, lsbc; the envelope triangle 60 trit, tril, ps, ts, rs, ftri, istore, estore and the fit cfenv[1..9]; 61 hstot, bl, teffl; the model hd, hl, hp, hpp, hr, hs, hs1, hs2, ht, htt, 62 hcomp[species, shell], lc; the Henyey work arrays hha, hhb, hhc, hsrf, 63 hmax, nhmax; tlumx[1..8], the luminosity by source; latmo, lderiv, 64 lmix, locond (for COEFFT); rotation fp, ft, eta2, omega, r0, 65 hg, hi, hjm, qiw, hkerot, hkerot0; itdone, the iterations done; 66 itlvl, the level of iteration; mxzone[12, 2], the convection zones; 67 cmixl, the mixing length. Arrays have index 0 unused and are changed 68 in place. 69 70 Returns (lfini, large, lnew, lreset, lsbc, ftri, istore, bl, teffl, 71 latmo, lderiv, lmix, locond, itdone): converged; corrections 72 too large; the flags and triangle data SURFBC updates; the new log L 73 and log Teff; the flags COEFFT sets; the iteration 74 count. 75 76 """ 77 78 global jcore, jenv, ddage, lflp 79 nml = nmlpar_mod 80 io = iounits_mod 81 82 # Fortran line 111. 83 if niter <= 0: 84 return (lfini, large, lnew, lreset, lsbc, ftri, istore, bl, teffl, 85 latmo, lderiv, lmix, locond, itdone) 86 87 # The luminosity, and zero counters (Fortran lines 112–117). 88 bl = math.log10(hl[m]) 89 kenv = 0 90 katm = 0 91 ksaha = 0 92 comp_mod.senv = hs[m] - hstot 93 94 # A new triangle for a new run or a reset (Fortran lines 118–121). 95 if lnew or (lreset and itlvl == 2): 96 lsbc = True 97 98 # Find new FP and FT if the model is rotating (Fortran lines 122–129). 99 if nml.lrot and lsbc: 100 fpl = fp[m] 101 ftl = ft[m] 102 else: 103 fpl = 1.0 104 ftl = 1.0 105 106 # Surface boundary condition: envelope composition, then SURFBC 107 # (Fortran lines 130–148). 108 if lsbc: 109 x = comp_mod.xnew 110 z = comp_mod.znew 111 plim = hp[m] 112 if nml.ldh: 113 debhu_mod.xxdh = hcomp[1, m] 114 debhu_mod.yydh = hcomp[2, m] + hcomp[4, m] 115 debhu_mod.zzdh = hcomp[3, m] 116 debhu_mod.zdh[1] = hcomp[5, m] + hcomp[6, m] 117 debhu_mod.zdh[2] = hcomp[7, m] + hcomp[8, m] 118 debhu_mod.zdh[3] = hcomp[9, m] + hcomp[10, m] + hcomp[11, m] 119 # SURFBC no longer sets TEFFL (since 7be4133). 120 (ftri, istore, lnew, lreset, ksaha, kenv, katm, 121 lflp) = surfbc(trit, tril, cfenv, ps, ts, rs, ftri, istore, estore, 122 lnew, lreset, ksaha, kenv, katm, hstot, bl, 123 teffl, x, z, fpl, ftl, lflp, plim, lc, m, cmixl) 124 125 # Reset the composition to the start of the step, then MIX (7/91) 126 # (Fortran lines 149–161). 127 if itlvl > 2 and delts > 0.0: 128 nspec = 11 129 if nml.lexcom: 130 nspec = 15 131 for i in range(1, m + 1): 132 for j in range(1, nspec + 1): 133 hcomp[j, i] = oldmod_mod.hcompp[j, i] 134 ddage, jcore, jenv = mix(delts, hcomp, hd, hl, hp, hr, hs, hs1, hs2, 135 hstot, ht, itlvl, lc, m, ddage, jcore, jenv, 136 mxzone, mxzon0, teffl, cmixl) 137 138 # Converged at the previous level and the triangle not flipped: done 139 # (Fortran lines 162–168; the call to OVROT is commented out). 140 if lsbc and not lflp and lfini: 141 return (lfini, large, lnew, lreset, lsbc, ftri, istore, bl, teffl, 142 latmo, lderiv, lmix, locond, itdone) 143 144 # The iterations (Fortran lines 169–310). 145 lfini = False 146 for iter_ in range(1, niter + 1): # loop 100 147 # Henyey reduction (Fortran lines 171–179). 148 hmax[1] = 0.0 149 hmax[2] = 0.0 150 hmax[3] = 0.0 151 hmax[4] = 0.0 152 # Suspected bug (see the page): jenv is unset or stale at levels 1 153 # and 2. 154 (latmo, lderiv, lmix, locond, 155 ksaha) = coefft(delts, m, hd, hha, hhb, hhc, hl, hmax, hp, hpp, hr, 156 hs, hs1, hs2, ht, htt, hcomp, lc, tlumx, latmo, 157 lderiv, lmix, locond, ksaha, fp, ft, 158 hkerot, hkerot0, jenv, teffl, cmixl) 159 160 # Renormalize TLUMX (not during the helium flash) (Fortran lines 161 # 180–189). 162 thl = (tlumx[1] + tlumx[2] + tlumx[3] + tlumx[4] + tlumx[5] 163 + tlumx[6] + tlumx[7] + tlumx[8]) 164 if not nml.lkuthe and thl > 0.0: 165 temp = hl[m] / thl 166 for j in range(1, 9): 167 tlumx[j] = tlumx[j] * temp 168 169 # Residuals of the P and T equations (not the boundary equations) 170 # (Fortran lines 190–202). 171 # Question (see the page): converged without solving. 172 if iter_ > 1: 173 if (hmax[1] <= nml.htoler[5, 1] and hmax[2] <= nml.htoler[5, 2] 174 and hmax[3] <= nml.htoler[5, 1] 175 and hmax[4] <= nml.htoler[5, 2]): 176 io.write(io.ishort, " R.H.S. BELOW TOLERANCES--P" 177 f"{hmax[1]:9.2E} T {hmax[2]:9.2E}" 178 f" R {hmax[3]:9.2E} L {hmax[4]:9.2E}") # format 20 179 io.write(io.ishort, f"{'':10}MODEL CONVERGED AFTER " 180 f"{itdone + 1:4d} ITERATIONS") # format 75 181 lfini = True 182 return (lfini, large, lnew, lreset, lsbc, ftri, istore, bl, 183 teffl, latmo, lderiv, lmix, locond, itdone) 184 185 # Henyey coefficients for the surface, from the envelope fit 186 # (Fortran lines 203–210). 187 hsrf[1] = -cfenv[1] 188 hsrf[2] = -cfenv[2] 189 hsrf[3] = cfenv[1] * hp[m] + cfenv[2] * ht[m] + cfenv[3] - hr[m] 190 temp = constants_mod.cln * hl[m] 191 hsrf[4] = -temp * cfenv[4] 192 hsrf[5] = -temp * cfenv[5] 193 hsrf[6] = temp * (cfenv[4] * hp[m] + cfenv[5] * ht[m] + cfenv[6] - bl) 194 195 # Back solution for the corrections (Fortran lines 211–212). 196 hsolve(m, hha, hhb, hsrf) 197 198 # Largest corrections, and where (Fortran lines 213–241). 199 for j in range(1, 5): 200 hmax[j] = abs(hhb[j, 1]) 201 nhmax[j] = 1 202 for i in range(2, m + 1): 203 test = abs(hhb[1, i]) 204 if hmax[1] <= test: 205 hmax[1] = test 206 nhmax[1] = i 207 test = abs(hhb[2, i]) 208 if hmax[2] <= test: 209 hmax[2] = test 210 nhmax[2] = i 211 test = abs(hhb[3, i]) 212 if hmax[3] <= test: 213 hmax[3] = test 214 nhmax[3] = i 215 # Question (see the page): the mixed absolute/relative L test. 216 test = min(abs(hhb[4, i]), abs(hhb[4, i] / hl[i])) 217 if hmax[4] <= test: 218 hmax[4] = test 219 nhmax[4] = i 220 hlmax = hmax[4] 221 222 # Converged below HTOLER(.,1); too large above HTOLER(.,2) (Fortran 223 # lines 242–252). 224 lfini = (hmax[1] < nml.htoler[1, 1] and hmax[2] < nml.htoler[2, 1] 225 and hmax[3] < nml.htoler[3, 1] 226 and hmax[4] < nml.htoler[4, 1]) 227 large = (hmax[1] > nml.htoler[1, 2] or hmax[2] > nml.htoler[2, 2] 228 or hmax[3] > nml.htoler[3, 2] or hmax[4] > nml.htoler[4, 2]) 229 for j in range(1, 5): 230 ij = nhmax[j] 231 hmax[j] = hhb[j, ij] 232 233 # The damping, raised each iteration; in the helium flash, by the 234 # size of the L correction (Fortran lines 253–259). 235 if nml.fcorr0 > 0.0: 236 ctol_mod.fcorr = min(1.0, ctol_mod.fcorr + nml.fcorri) 237 fact = ctol_mod.fcorr 238 if nml.lkuthe: 239 if hlmax <= 5.0e-1: 240 fact = 8.0e-1 241 if hlmax <= 5.0e-3: 242 fact = 1.0 243 244 # The corrections printed; stop if too large (Fortran lines 245 # 260–280). 246 tlumh = tlumx[1] + tlumx[2] + tlumx[3] + tlumx[4] 247 tlumhe = tlumx[5] + tlumx[8] 248 if nml.lcorr: 249 io.write(io.ishort, f" DEL-L/L {_l(lfini):>2}{hmax[4]:12.4E}" 250 f"{nml.htoler[4, 1]:12.4E}{'':5}{nhmax[4]:5d}") # fmt 60 251 io.write(io.ishort, " CORR" 252 f"{nhmax[1]:5d}P{hmax[1]:9.2E}{nhmax[2]:5d}T{hmax[2]:9.2E}" 253 f"{nhmax[3]:5d}R{hmax[3]:9.2E}{nhmax[4]:5d}L{hmax[4]:9.2E}" 254 f" F={fact:5.3f} E-HY{tlumh:10.3E} HE{tlumhe:10.3E}" 255 f" NU{tlumx[6]:10.3E} G{tlumx[7]:10.3E}") # format 70 256 if lfini: 257 io.write(io.ishort, f"{'':10}MODEL CONVERGED AFTER " 258 f"{itdone + 1:4d} ITERATIONS") # format 75 259 if large: 260 io.write(io.iowr, " -----CORRECTIONS EXCEEDED TOLERANCES") 261 io.write(io.ishort, " -----CORRECTIONS EXCEEDED TOLERANCES") 262 return (lfini, large, lnew, lreset, lsbc, ftri, istore, bl, 263 teffl, latmo, lderiv, lmix, locond, itdone) 264 265 # Apply the corrections, damped by fact (Fortran lines 281–291). 266 for i in range(1, m + 1): 267 temp = fact * hhb[1, i] 268 hp[i] = hp[i] + temp 269 hpp[i] = hpp[i] + temp 270 temp = fact * hhb[2, i] 271 ht[i] = ht[i] + temp 272 htt[i] = htt[i] + temp 273 hr[i] = hr[i] + fact * hhb[3, i] 274 hl[i] = hl[i] + fact * hhb[4, i] 275 276 # The new log L and log Teff, and the rotation (Fortran lines 277 # 292–307). 278 bl = math.log10(hl[m]) 279 teffl = cfenv[7] * hp[m] + cfenv[8] * ht[m] + cfenv[9] 280 if nml.lrot: 281 # The convective flags LCZO from OVROT, stored by HPOINT (or 282 # STARIN, for the starting model) (MHP via Claude, 10/29/26). 283 getrot(hd, hjm, hr, hs, hs2, oldmod_mod.lczo, m, eta2, hi, omega, 284 qiw, r0) 285 fpft(hd, hr, hs, m, omega, eta2, fp, ft, hg, r0) 286 for i in range(1, m + 1): 287 hj = hjm[i] * hs2[i] 288 hkerot[i] = 0.5 * omega[i] * hj 289 290 # Fortran lines 308–310. 291 itdone = itdone + 1 292 if lfini: 293 break 294 return (lfini, large, lnew, lreset, lsbc, ftri, istore, bl, teffl, 295 latmo, lderiv, lmix, locond, itdone)