bursmix.f

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

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

Variables on this page (11)

hcompp, json, lexcom, pscal0, pscal1, rbe9, rbe90, rli6, rli60, rli7, rli70

Purpose

This file evolves the abundances over one time step when rotational mixing and nuclear burning act together. In its own words, it alternates "between nuclear burning plus gravitational settling and mixing in a series of progressively smaller timesteps, extrapolating to the answer for zero timestep." The diffusion coefficients, structure and reaction rates are held fixed during the step.

How it fits

Inputs: the arguments, the model (HD, HL, HP, HR, HS, HS1, HS2, HT, HV), the mixing coefficients COD2, the time step DELTS, the zones (MRZONE, MXZONE, NRZONE, NZONE, LCZ); LEXCOM from NMLPAR_MOD (whether the extended set of 15 species, with the light elements, is followed); HCOMPP in OLDMOD_MOD, the composition at the start of the step.
Who calls it: GETW, the rotation step, only when LINSTB is true, LALLCZ false and LBURS true (getw.f, lines 455–460 ↗). GETW sets its local LBURS to false on every call and never changes it (getw.f, line 140 ↗), so at 7be4133 BURSMIX is never called, nor BSROTMIX through it: unused?
What it calls: ROTMIX (burning, by its header with KEMCOM in radiative regions), NDIFCOM (mixing by angular-momentum transport), LIBURN2 (light-element burning, if LEXCOM), and BSROTMIX (the extrapolation).
Outputs: HCOMP, the abundances at the end of the step. If LEXCOM, it copies the new lithium and beryllium rates (NEWRAT_MOD) into the old ones (OLDRAT_MOD), and sets JENV0 = JENV1, PSCAL0 = PSCAL1 (LIOV_MOD).
Modules: OLDMOD_MOD (18 units, among them CRRECT, KEMCOM, LIBURN); NEWRAT_MOD, OLDRAT_MOD, LIOV_MOD (GETW, LIBURN, LIBURN2, LIRATE88); NMLPAR_MOD.

The method

Save the composition at the start of the step (HCOMPP), and $X$, $Y$, $Z$ of the input HCOMP, which the caller's full step has already produced (lines 75–82 ↗).
Pass the composition of the full step, already computed by the caller, to BSROTMIX as sequence 1.
For sequences $k = 2, \dots, 11$, with $n_k$ = 2, 4, 6, …, 20 substeps of $\Delta t/n_k$: reset the composition, then for each substep burn and, if LDIFY, settle (ROTMIX, which calls MICRODIFF or GRSETT), mix (NDIFCOM) and, if LEXCOM, burn the light elements (LIBURN2); the end of each substep becomes the start of the next (HCOMPP, which KEMCOM reads as the start-of-step state). Then call BSROTMIX; stop when it reports convergence (lines 87–126 ↗). On convergence HCOMPP is set to the end-of-step composition (lines 118–123 ↗), and otherwise it is left there by the last substep, though the comment says it is overwritten "temporarily" (lines 72–74 ↗); GETW does the same on its other path, copying the end-of-step composition into HCOMPP when LBURS is false (getw.f, lines 372–377 ↗), so it is the start of the next step (settled from the code, Nov 9, 2026).

Question 1 (line 52 ↗; lines 95–96 ↗). Two NSEQ sequences. NSEQ here starts 1, 2, 4, 6, … (step by 2), while BSSTEP uses 2, 4, 6, 8, 12, … Is the choice deliberate?
If no sequence converges, print a warning and keep the composition of the last (finest) sequence.
Restore $X$ and $Z$ to their input values (those of the caller's full step, not those at the start of the time step), and set $Y = 1 - X - Z - X_3$ (lines 131–136 ↗).

Question 2 (lines 132–136 ↗). $X$ and $Z$ restored. After the extrapolation, $X$ and $Z$ are reset to their values on entry (HKEEP), that is, to the caller's full step, so the burning and mixing of hydrogen computed in the substeps are discarded, and $Y$ absorbs the difference. Is this reset meant?

Only species 4 onward (from $^3$He) are extrapolated (JBEG = 4).

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. (cannot be reached) NDIFCOM gets a region never set when `LOK` is false.

(lines 28–31 ↗; lines 46–49 ↗; lines 103–106 ↗) IBEG, IEND, ECOD2 and EM are locals of BURSMIX, listed as input variables but set nowhere in it; only NDIFCOM sets them, through MIXGRID, and only when LOK is true. With LOK false, NDIFCOM passes them straight to MIXCOM (ndifcom.f, lines 64–66 ↗). GETW passes the LOK that SECULR returns, which is false when SECULR finds no instability and leaves early (seculr.f, line 202 ↗, seculr.f, line 289 ↗). What follows: On the first such call IBEG = IEND = 0 if the compiler zeroes saved locals, and MIXCOM reads HCOMP(J,0), out of bounds (mixcom.f, line 61 ↗), before finding nothing to do; after an earlier call with LOK true, MIXCOM mixes again the last region of that call, with its old coefficients, in a model SECULR found stable.

A possible fix:

bursmix.f

@@ lines 100–109 @@
             CALL ROTMIX(DT,HCOMP,HS2,HT,M,MRZONE,MXZONE,
      *               NRZONE,NZONE,HSTOT,HD,HS,HR,HP,LCZ,HS1)
 C PERFORM MIXING
-            CALL NDIFCOM(DT,COD2,ECOD2,EM,HD,HL,HP,HR,HS,HS1,HS2,
+C 11/09/26 KC via Claude
+C code changed to mix only when SECULR found unstable regions (LOK);
+C IBEG, IEND, ECOD2 and EM are set only by such a call.
+            IF(LOK)THEN
+               CALL NDIFCOM(DT,COD2,ECOD2,EM,HD,HL,HP,HR,HS,HS1,HS2,
 C      *                   HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M,NTOT,  ! KC 2025-05-31
-     *                   HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M,
-     *                   HCOMP,JBEG,JEND)
+     *                      HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M,
+     *                      HCOMP,JBEG,JEND)
+            ENDIF
             IF(LEXCOM)THEN
 C                CALL LIBURN2(DT,HCOMP,HD,HR,HS1,HS2,HT,JENV1,JENV0,M)  ! KC 2025-05-31
                CALL LIBURN2(DT,HCOMP,HR,HS1,HS2,HT,JENV1,JENV0,M)

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Cannot be reached at 7be4133. GETW sets LBURS = .FALSE. on every call (getw.f, line 140 ↗) and calls BURSMIX only when it is true (getw.f, lines 455–461 ↗); no namelist can change that. No counting line fired in any case run, the five rotating test cases run (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) among them. A run would reach it only if the extrapolation were switched on in GETW. Intent: none found.

Definitions

subroutineBURSMIX

lines 2–147 · GitHub

Burn–mix substepping over one time step for increasing numbers of substeps, extrapolated to zero substep by BSROTMIX; resets $X$, $Z$ and updates $Y$.

Uses modules DIMENSIONS_MOD, LIOV_MOD, NEWRAT_MOD, NMLPAR_MOD, OLDMOD_MOD, OLDRAT_MOD

Called by GETW

Flowchart of its calls, in order (5)

BSROTMIX 86

for III = 2, NBURS, lines 87–126

for JJJ = 1,NSTEP, lines 97–116

ROTMIX 100→NDIFCOM 103

if (LEXCOM), lines 107–110

LIBURN2 109

BSROTMIX 117

Change log

The changes since rp_1

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

1C $$$$$$ 2 SUBROUTINE BURSMIX(COD2,DELTS,HCOMP,HD,HL,HP,HR,HS,HS1,HS2,HSTOT, 3 * * HT,HV,IMAX,IMIN,JENV0,JENV1,LOK,LCZ,M,MRZONE, 4 * * MXZONE,NRZONE,NZONE) 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/OLDRAT/ into OLDRAT_MOD. 11 USE OLDRAT_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/OLDMOD/ into OLDMOD_MOD. 14 USE OLDMOD_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/NEWRAT/ into NEWRAT_MOD. 17 USE NEWRAT_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/LIOV/ into LIOV_MOD. 20 USE LIOV_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 23 USE NMLPAR_MOD 24C 10/20/26 MHP via Claude 25C code changed to replace implicit typing with IMPLICIT NONE and 26C explicit declarations. 27 IMPLICIT NONE 28 REAL*8 COD2,DELTS,DT,ECOD2,EM,HCOMP,HCSAV,HD,HKEEP,HL,HP,HR,HS, 29 * * HS1,HS2,HSTOT,HT,HV 30 INTEGER I,IBEG,IEND,III,IMAX,IMIN,J,JBEG,JEND,JENV0,JENV1,JJJ,M, 31 * * MRZONE,MXZONE,NBURS,NRZONE,NSEQ,NSPEC,NSTEP,NZONE 32 LOGICAL*4 LCONV,LCZ,LOK 33C 10/12/26 MHP via Claude 34C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 35C namelist variables that the physics routines read from 27 COMMON 36C blocks which also held derived or working members. Those blocks keep 37C only their non-namelist members; five of them -- NEWMX, VMULT, 38C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 39C Derivation-only namelist variables from the same blocks went to 40C COMMON/NMLBYPASS/ instead, not here. 41C 42C Declared identically in all 121 files that hold any of these members. 43C Generated from one member list, so every copy is byte-identical; the 44C linker's 'size of symbol changed' warning is the check that they are. 45C INPUT VARIABLES 46 DIMENSION COD2(JSON),HCOMP(15,JSON),HD(JSON),HL(JSON), 47 * * HP(JSON),HR(JSON),HS(JSON),HS1(JSON),HS2(JSON),HT(JSON), 48 * * HV(JSON),LCZ(JSON),MRZONE(13,2),MXZONE(12,2), 49 * * ECOD2(JSON),EM(JSON) 50C LOCAL VARIABLES 51 DIMENSION HCSAV(15,JSON),NSEQ(11),HKEEP(3,JSON) 52 DATA NSEQ/1,2,4,6,8,10,12,14,16,18,20/ 53 SAVE 54C THIS SUBROUTINE IS DESIGNED TO ALTERNATE BETWEEN NUCLEAR BURNING 55C PLUS GRAVITATIONAL SETTLING AND MIXING IN A SERIES OF PROGRESSIVELY 56C SMALLER TIMESTEPS, EXTRAPOLATING TO THE ANSWER FOR ZERO TIMESTEP. 57C IT ASSUMES THAT THE DIFFUSION COEFFICIENTS, STRUCTURE, AND NUCLEAR REACTION 58C RATES ARE HELD CONSTANT DURING THE TIME STEP. THE PROCEDURE IS AS 59C FOLLOWS: 60C STORE THE RESULTS FOR THE WHOLE TIMESTEP IN ENTRY ONE. 61C THEN CALL THE NUCLEAR BURNING ROUTINES TO BURN AND THE ROTATIONAL 62C MIXING ROUTINE TO MIX FOR A SERIES OF SMALLER TIME STEPS, STORING THE 63C ANSWER AS A FUNCTION OF TIME STEP. 64 IF(LEXCOM)THEN 65 NSPEC = 15 66 ELSE 67 NSPEC = 11 68 ENDIF 69 JBEG = 4 70 JEND = NSPEC 71 NBURS = 11 72C STORE ALL ELEMENTS HERE - HCOMPP IS USED AS THE START OF 73C TIMESTEP VARIABLE IN KEMCOM, SO IT NEEDS TO BE TEMPORARILY 74C OVERWRITTEN IN THE SOLUTION PROCESS 75 DO I = 1,M 76 DO J = 1,NSPEC 77 HCSAV(J,I) = HCOMPP(J,I) 78 END DO 79 DO J = 1,3 80 HKEEP(J,I) = HCOMP(J,I) 81 END DO 82 END DO 83C STORE RESULTS OF INITIAL CALCULATION AS FIRST ENTRY IN BS 84C ROUTINE 85 III = 1 86 CALL BSROTMIX(DELTS,HCOMP,III,M,JBEG,NSPEC,NSEQ,LCONV) 87 DO III = 2, NBURS 88C RESET COMPOSITION TO INITIAL VALUE. 89 DO I = 1,M 90 DO J = 1,NSPEC 91 HCOMP(J,I) = HCSAV(J,I) 92 HCOMPP(J,I) = HCSAV(J,I) 93 END DO 94 END DO 95 NSTEP = NSEQ(III) 96 DT = DELTS/DFLOAT(NSTEP) 97 DO JJJ = 1,NSTEP 98C PERFORM NUCLEAR BURNING 99C CALL ROTMIX(DT,HCOMP,HS2,HT,5,M,MRZONE,MXZONE, ! KC 2025-05-31 100 CALL ROTMIX(DT,HCOMP,HS2,HT,M,MRZONE,MXZONE, 101 * * NRZONE,NZONE,HSTOT,HD,HS,HR,HP,LCZ,HS1) 102C PERFORM MIXING 103 CALL NDIFCOM(DT,COD2,ECOD2,EM,HD,HL,HP,HR,HS,HS1,HS2, 104C * HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M,NTOT, ! KC 2025-05-31 105 * * HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M, 106 * * HCOMP,JBEG,JEND) 107 IF(LEXCOM)THEN 108C CALL LIBURN2(DT,HCOMP,HD,HR,HS1,HS2,HT,JENV1,JENV0,M) ! KC 2025-05-31 109 CALL LIBURN2(DT,HCOMP,HR,HS1,HS2,HT,JENV1,JENV0,M) 110 ENDIF 111 DO I = 1,M 112 DO J = 1,NSPEC 113 HCOMPP(J,I) = HCOMP(J,I) 114 END DO 115 END DO 116 END DO 117 CALL BSROTMIX(DELTS,HCOMP,III,M,JBEG,NSPEC,NSEQ,LCONV) 118 IF(LCONV)THEN 119 DO I = 1,M 120 DO J = 1,NSPEC 121 HCOMPP(J,I) = HCOMP(J,I) 122 END DO 123 END DO 124 GOTO 10 125 ENDIF 126 END DO 127C FAILED TO CONVERGE; PRINT WARNING 128C IN THIS CASE THE UNEXTRAPOLATED FINAL COMPOSITION IS USED. 129 WRITE(*,5) 130 5 5 FORMAT(' WARNING - EXTRAPOLATION IN BSBURN DID NOT CONVERGE') 131 10 10 CONTINUE 132 DO I = 1,M 133 HCOMP(1,I) = HKEEP(1,I) 134 HCOMP(3,I) = HKEEP(3,I) 135 HCOMP(2,I) = 1.0D0 - HCOMP(1,I) - HCOMP(3,I) - HCOMP(4,I) 136 END DO 137 IF(LEXCOM)THEN 138 DO I = 1,M 139 RLI60(I) = RLI6(I) 140 RLI70(I) = RLI7(I) 141 RBE90(I) = RBE9(I) 142 END DO 143 JENV0 = JENV1 144 PSCAL0 = PSCAL1 145 ENDIF 146 RETURN 147 END

1"""Python version of bursmix.f (YREC 7be4133): burning and mixing over one 2time step in progressively smaller substeps, extrapolated to zero substep.""" 3 4import numpy as np 5 6from pyrec import dimensions_mod 7from pyrec import liov_mod 8from pyrec import newrat_mod 9from pyrec import nmlpar_mod 10from pyrec import oldmod_mod 11from pyrec import oldrat_mod 12from pyrec.bsrotmix import bsrotmix 13from pyrec.liburn2 import liburn2 14from pyrec.ndifcom import ndifcom 15from pyrec.rotmix import rotmix 16 17JSON = dimensions_mod.json 18 19# Locals kept between calls (SAVE). ECOD2, EM, IBEG and IEND are filled 20# only by NDIFCOM (through MIXGRID) when LOK is true. 21# Suspected bug (see the page): with LOK false they are read as left here. 22ecod2 = np.zeros(JSON + 1) 23em = np.zeros(JSON + 1) 24ibeg = 0 25iend = 0 26hcsav = np.zeros((16, JSON + 1)) 27hkeep = np.zeros((4, JSON + 1)) 28# DATA NSEQ: the substep numbers of the 11 sequences (index 0 unused). 29# Question (see the page): 1, 2, 4, 6, ..., unlike bsstep's 2, 4, 6, 8, 12, ... 30nseq = np.array([0, 1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]) 31 32 33def bursmix(cod2: np.ndarray, delts: float, hcomp: np.ndarray, hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, 34 hr: np.ndarray, hs: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, hstot: float, ht: np.ndarray, 35 hv: np.ndarray, imax: int, imin: int, jenv0: int, jenv1: int, lok: bool, lcz: np.ndarray, m: int, 36 mrzone: np.ndarray, mxzone: np.ndarray, nrzone: int, nzone: int) -> int: 37 """Alternate nuclear burning (with settling) and rotational mixing in 38 sequences of 2, 4, ..., 20 substeps, and extrapolate the abundances 39 from species 4 on to zero substep with bsrotmix; then restore X and Z 40 to their values on entry and set Y = 1 - X - Z - X(He3). 41 42 Translates: bursmix.f#BURSMIX (7be4133) 43 44 Arguments: the mixing coefficients cod2; the time step delts (s); the 45 composition hcomp[species, shell] after the caller's full step 46 (changed in place); the model hd, hl, hp, hr, hs, hs1, hs2, hstot, 47 ht, hv; imax, imin; jenv0, jenv1 (for the light elements); lok (from 48 SECULR); lcz; m shells; the zones mrzone, mxzone, nrzone, nzone. 49 oldmod_mod.hcompp is overwritten and left at the end-of-step 50 composition. 51 52 Returns jenv0: set to jenv1 when nmlpar_mod.lexcom, else as passed in. 53 54 """ 55 56 global ibeg, iend 57 n = nmlpar_mod 58 hcompp = oldmod_mod.hcompp 59 60 # The species extrapolated, from He3 on (Fortran lines 64–71). 61 if n.lexcom: 62 nspec = 15 63 else: 64 nspec = 11 65 jbeg = 4 66 jend = nspec 67 nburs = 11 68 69 # Store all elements: HCOMPP is the start-of-step composition in 70 # KEMCOM, so it is overwritten during the substeps (Fortran lines 75–82). 71 for i in range(1, m + 1): 72 for j in range(1, nspec + 1): 73 hcsav[j, i] = hcompp[j, i] 74 for j in range(1, 4): 75 hkeep[j, i] = hcomp[j, i] 76 77 # The caller's full step is the first entry of the extrapolation (Fortran lines 85–86). 78 iii = 1 79 lconv = False 80 lconv = bsrotmix(delts, hcomp, iii, m, jbeg, nspec, nseq, lconv) 81 converged = False 82 83 # The sequences of substeps (Fortran lines 87–126). 84 for iii in range(2, nburs + 1): 85 # Reset the composition to its initial value. 86 for i in range(1, m + 1): 87 for j in range(1, nspec + 1): 88 hcomp[j, i] = hcsav[j, i] 89 hcompp[j, i] = hcsav[j, i] 90 nstep = int(nseq[iii]) 91 dt = delts / float(nstep) 92 for jjj in range(1, nstep + 1): 93 # Nuclear burning (and settling). 94 rotmix(dt, hcomp, hs2, ht, m, mrzone, mxzone, nrzone, nzone, 95 hstot, hd, hs, hr, hp, lcz, hs1) 96 # Mixing. Suspected bug (see the page): with lok False, ibeg, iend, 97 # ecod2 and em were never set. 98 ibeg, iend = ndifcom(dt, cod2, ecod2, em, hd, hl, hp, hr, hs, 99 hs1, hs2, hstot, hv, ibeg, iend, imax, 100 imin, lcz, lok, m, hcomp, jbeg, jend) 101 if n.lexcom: 102 liburn2(dt, hcomp, hr, hs1, hs2, ht, jenv1, jenv0, m) 103 for i in range(1, m + 1): 104 for j in range(1, nspec + 1): 105 hcompp[j, i] = hcomp[j, i] 106 lconv = bsrotmix(delts, hcomp, iii, m, jbeg, nspec, nseq, lconv) 107 if lconv: 108 for i in range(1, m + 1): 109 for j in range(1, nspec + 1): 110 hcompp[j, i] = hcomp[j, i] 111 converged = True 112 break # GOTO 10 113 114 # Failed to converge: the unextrapolated final composition is used (Fortran lines 129–130). 115 if not converged: 116 print(' WARNING - EXTRAPOLATION IN BSBURN DID NOT CONVERGE') 117 118 # label 10. X and Z restored, Y by closure (Fortran lines 131–136). 119 # Question (see the page): X and Z restored to their values on entry. 120 for i in range(1, m + 1): 121 hcomp[1, i] = hkeep[1, i] 122 hcomp[3, i] = hkeep[3, i] 123 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 124 125 # The light-element rates and envelope base for the next step (Fortran lines 137–145). 126 if n.lexcom: 127 for i in range(1, m + 1): 128 oldrat_mod.rli60[i] = newrat_mod.rli6[i] 129 oldrat_mod.rli70[i] = newrat_mod.rli7[i] 130 oldrat_mod.rbe90[i] = newrat_mod.rbe9[i] 131 jenv0 = jenv1 132 liov_mod.pscal0 = liov_mod.pscal1 133 return jenv0

bursmix.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (cannot be reached) NDIFCOM gets a region never set when LOK is false.

Definitions

subroutineBURSMIX

Change log

1. (cannot be reached) NDIFCOM gets a region never set when `LOK` is false.