dburn.f

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

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

Variables on this page (8)

compacc, dmdt0, drate, drate0, fmassacc, hcompp, json, lmdot

Purpose

This file burns deuterium over a time step in one shell or one mixed region (a convection zone), with the rate varying linearly from its start-of-step to its end-of-step value, and updates hydrogen and $^3$He accordingly. With mass accretion, it also burns the accreted deuterium in the surface convection zone.

How it fits

Inputs: the arguments IBEGIN, IEND (the shell or zone), M, HS2 (the mass of each shell), HCOMP, DT (the time step; MIX passes $\Delta t$ in Gyr, DDAGE*1.0D-9); the start-of-step composition HCOMPP (OLDMOD_MOD); the rates DRATE0, DRATE (DEUTER_MOD, from DEUTRATE) and FMASSACC (DEUTER_MOD, the mass accreted over the step in grams, set by MDOT when mass is gained); the settings LMDOT, DMDT0, COMPACC (NMLPAR_MOD, formerly /MASSCHG/).
Who calls it: MIX, only with LEXCOM (the extended composition, 15 species): for each radiative shell, from the centre out until $T$ falls below TCUT(1), and for every convection zone.
Output: HCOMP (species 1, 4 and 12: H, $^3$He, D) over the shell or zone.

The method

For a zone, mass-weighted means of $X$, $X_D$, $X_3$ and of the rates, with the zone's mass SUMM; for a single shell, SUMM is the shell's own mass.
If $X_D$ (with the accreted D, at the surface) is below $10^{-11}$, set it to zero and return (leaving $X$ and $^3$He as they are).
Burn in ISTEP = $\lfloor r\Delta t/0.1\rfloor + 1$ substeps, the rate rising linearly, with $r_k$ the rate at the substep's midpoint: $$X_D \leftarrow X_D\,e^{-2 r_k\,\Delta t/n}.$$ The code does not say what the factor 2 is for; DRATE is divided by the deuteron's mass AMDEUT ("Avagadros number divided by the mass of the deuteron in amu," DEUTRATE), so the second question below asks.

Question 1 (lines 89–94 ↗). The substep size. The comment says the burning per substep is limited to 0.5 in the natural log; the code limits $r\,\Delta t/n$ to 0.1, that is $2r\Delta t/n \le 0.2$. Is only the comment off, or was a coarser substep (0.25 for $r\Delta t/n$) meant?

Question 2 (line 102 ↗; lines 120–122 ↗; deutrate.f, line 62 ↗; deutrate.f, line 81 ↗). The factor 2 in the burning. The substeps use $X_D \leftarrow X_D e^{-2r\Delta t}$. DRATE carries $N_A/A_D$ with $A_D \approx 2$, so $2\,$DRATE $\approx X\rho N_A\langle\sigma v\rangle$ per Gyr, which is the destruction rate per deuteron by protons: is that what the factor 2 is for? If so, $N_A/A_p$ with $A_p \approx 1$ would be the natural factor, and the 2 and AMDEUT differ slightly ($2/2.013553$); DBURNM uses the same factor. 4. Mass bookkeeping of D + p $\to$ $^3$He: with $\Delta X_D < 0$ the change of deuterium, $$\Delta X = \tfrac12\,\Delta X_D, \qquad \Delta X_3 = -\tfrac32\,\Delta X_D.$$
Accretion (LMDOT, surface zone): matter accreted earlier in the step burns fully, new matter is exposed to half a substep's burning; then the zone's values are mass-weighted means of the old and the accreted matter.

FMASSACC is a mass in grams, not a fraction. The comment defines it as "DMDT*DT/original CZ mass" (lines 107–109 ↗), but MDOT sets it to DELM = DMDT*CMSUN*DELTS/CSECYR, the mass accreted over the step in grams (mdot.f, line 159 ↗, mdot.f, line 320 ↗). DBURN's means, weighted by SUMM (a sum of HS2, in grams: lines 48–62 ↗) and FMASSACC, are therefore consistent (and DBURNM's are since a8a5950); the same stale comment stands in DBURNM (dburnm.f, lines 111–112 ↗) and LIBURN (liburn.f, line 451 ↗).

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. (harmless in the runs) `FMASSACC` is read under mass loss.

(line 75 ↗; line 106 ↗; lines 128–132 ↗) The accretion branch is taken whenever LMDOT is set and IEND = M, also with mass loss ($\dot M \le 0$), while the zeroing test above it checks DMDT0 > 0 too. MDOT sets FMASSACC only when mass is gained (mdot.f, line 300 ↗, mdot.f, line 320 ↗), and DEUTER_MOD gives it no initial value (deuter_mod.f, line 30 ↗). In a run that only loses mass, the blend therefore uses a never-set FMASSACC: Compilers usually start module variables at zero, and with zero the branch gives the same result as the ELSE branch, so the code is correct in practice but rests on an unset variable. Only with LMDOT (off by default).

A possible fix:

dburn.f

@@ lines 103–109 @@
          R0 = R0 + DELRATE
       END DO
 C INCLUDE MASS ACCRETION FROM DEUTERIUM BURNING
-      IF(LMDOT .AND. IEND.EQ.M)THEN
+C 11/09/26 KC via Claude
+C code changed to take the accretion branch only when mass is gained,
+C as in the zeroing test above, since MDOT sets FMASSACC only then.
+      IF(LMDOT .AND. IEND.EQ.M .AND. DMDT0.GT.0.0D0)THEN
 C ACCRETED MATTER IS EXPOSED TO BURNING FOR, ON
 C AVERAGE. 1/2 OF THE TIMESTEP.  FMASSACC IS
 C DEFINED AS DMDT*DT/ORIGINAL CZ MASS.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): 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), only a made-up mass-loss case took the branch: YREC's testsuite/Test_m3p0feh+0p0GS98_std25_TAMS with LMDOT and LREIMER on. There it ran 8,341 times, every time with DMDT0 = $-10^{-14}$ (mass lost), so FMASSACC was read unset; it was 0 each time. With the fix, the track (8,363 models, as unfixed; the run stopped after the mass fell to 0.11 $M_\odot$) moves at its end by up to $6.2\times10^{-7}$ dex in $\log L$ and $4.1\times10^{-7}$ in $\log T_c$, against $2\times10^{-7}$ and $3\times10^{-8}$ for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$): $\log L$ within the noise, $\log T_c$ above it but tiny, probably rounding (not traced), since with FMASSACC = 0 the two branches compute the same values by different operations. The other made-up mass-loss case (YREC's 1 $M_\odot$ case to helium ignition with LMDOT and LREIMER on) did not take it. As read, the zero start value makes the result right in practice; fragile code.

Definitions

subroutineDBURN

lines 4–153 · GitHub

Deuterium burning over the step in a shell or zone, with a linearly varying rate in substeps, the H and $^3$He bookkeeping, and accreted deuterium at the surface.

Uses modules DEUTER_MOD, DIMENSIONS_MOD, NMLPAR_MOD, OLDMOD_MOD

Called by MIX ×2

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE, and the unused /CONST1/ gone.

Changed since a8a5950. For a single shell, the region mass SUMM is now set to the shell's own mass; before, the accretion blend used the value left from an earlier call (SAVE). This settles former question 2 and, with the fix of DBURNM's weighting, the second part of question 1; the later questions are renumbered.

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C DBURN - COMPUTE THE ABUNDANCE CHANGES RESULTING FROM DEUTERIUM BURNING 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4 SUBROUTINE DBURN(IBEGIN,IEND,M,HS2,HCOMP,DT) 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/OLDMOD/ into OLDMOD_MOD. 11 USE OLDMOD_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/DEUTER/ into DEUTER_MOD. 14 USE DEUTER_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 17 USE NMLPAR_MOD 18C 10/20/26 MHP via Claude 19C code changed to replace implicit typing with IMPLICIT NONE and 20C explicit declarations. 21 IMPLICIT NONE 22 REAL*8 DELRATE,DT,DXH2,DXH2A,DXH2O,FSTEP,HCOMP,HS2,R0,RATE,RATE0, 23 * * RATE1,SUMM,SUMR0,SUMR1,TSTEP,X,X2NEW,XH2,XH2N,XH2NA, 24 * * XH2NACC,XH2NO,XH2TEST,XHE3,XHE3N,XN 25 INTEGER I,IBEGIN,IEND,ISTEP,M 26C 10/13/26 MHP via Claude: removed the unused COMMON/CONSTANTS/ declaration 27C (no member referenced here; residue of the obsolete-option removals). 28C 10/12/26 MHP via Claude 29C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 30C namelist variables that the physics routines read from 27 COMMON 31C blocks which also held derived or working members. Those blocks keep 32C only their non-namelist members; five of them -- NEWMX, VMULT, 33C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 34C Derivation-only namelist variables from the same blocks went to 35C COMMON/NMLBYPASS/ instead, not here. 36C 37C Declared identically in all 121 files that hold any of these members. 38C Generated from one member list, so every copy is byte-identical; the 39C linker's 'size of symbol changed' warning is the check that they are. 40C DATA C21/5.240358E-8/ 41 DIMENSION HS2(JSON),HCOMP(15,JSON) 42 SAVE 43 IF(IBEGIN.EQ.IEND)THEN 44C 10/29/26 MHP via Claude 45C A single shell's region mass is its own mass. SUMM was set only in 46C the multi-shell branch, so the accretion blend below used the previous 47C call's region mass (SAVEd) for a single shell. 48 SUMM = HS2(IBEGIN) 49 X = HCOMPP(1,IBEGIN) 50 XH2 = HCOMPP(12,IBEGIN) 51 XHE3 = HCOMPP(4,IBEGIN) 52 RATE0 = DRATE0(IBEGIN) 53 RATE1 = DRATE(IBEGIN) 54 ELSE 55 SUMM = 0.0D0 56 SUMR0 = 0.0D0 57 SUMR1 = 0.0D0 58 X = 0.0D0 59 XH2 = 0.0D0 60 XHE3 = 0.0D0 61 DO I = IBEGIN,IEND 62 SUMM = SUMM + HS2(I) 63 SUMR0 = SUMR0 + HS2(I)*DRATE0(I) 64 SUMR1 = SUMR1 + HS2(I)*DRATE(I) 65 X = X + HCOMPP(1,I)*HS2(I) 66 XH2 = XH2 + HCOMPP(12,I)*HS2(I) 67 XHE3 = XHE3 + HCOMPP(4,I)*HS2(I) 68 END DO 69 RATE0 = SUMR0/SUMM 70 RATE1 = SUMR1/SUMM 71 X = X/SUMM 72 XH2 = XH2/SUMM 73 XHE3 = XHE3/SUMM 74 ENDIF 75 IF(LMDOT .AND. IEND.EQ.M .AND. DMDT0.GT.0.0D0)THEN 76 XH2TEST = (XH2*SUMM+COMPACC(12)*FMASSACC)/(SUMM+FMASSACC) 77 ELSE 78 XH2TEST = XH2 79 ENDIF 80 IF(XH2TEST.LT.1.0D-11)THEN 81 DO I = IBEGIN,IEND 82 HCOMP(12,I) = 0.0D0 83 END DO 84 GOTO 9999 85 ENDIF 86C BURN DEUTERIUM IN A SERIES OF SMALL STEPS, 87C INCREMENTING THE BURNING RATE FROM THE STARTING 88C TO ENDING VALUES. 89C RESTRICT THE MAXIMUM BURNING TO 0.5 IN THE 90C NATURAL LOG (<1% ERROR FROM A SIMPLE SUM 91C OF DISCRETE STEPS OF THIS SIZE). 92C MHP 10/02 INITIALIZE RATE 93 RATE = MAX(RATE0,RATE1) 94 ISTEP = INT(RATE*DT/0.1D0)+1 95 FSTEP = DFLOAT(ISTEP) 96 TSTEP = DT/FSTEP 97 R0 = RATE0 98 DELRATE = (RATE1-RATE0)/FSTEP 99 XH2NO = XH2 100 DO I = 1,ISTEP 101 RATE = R0 + 0.5D0*DELRATE 102 XH2NO = XH2NO*EXP(-2.0D0*RATE*TSTEP) 103 R0 = R0 + DELRATE 104 END DO 105C INCLUDE MASS ACCRETION FROM DEUTERIUM BURNING 106 IF(LMDOT .AND. IEND.EQ.M)THEN 107C ACCRETED MATTER IS EXPOSED TO BURNING FOR, ON 108C AVERAGE. 1/2 OF THE TIMESTEP. FMASSACC IS 109C DEFINED AS DMDT*DT/ORIGINAL CZ MASS. 110C BURN BOTH THE ACCRETED D AND THE ORIGINAL D 111C SEPARATELY AND FIND THE NEW MASS-WEIGHTED 112C AVERAGE D. NOTE THAT ACCRETION OF ELEMENTS 113C 1-11 HAS ALREADY BEEN TREATED. 114 R0 = RATE0 115 XH2NACC = COMPACC(12) 116 XH2NA = 0.0D0 117 DO I = 1,ISTEP 118 RATE = R0 + 0.5D0*DELRATE 119C FULLY BURN PRIOR ACCRETED MATTER 120 XH2NA = XH2NA*EXP(-2.0D0*RATE*TSTEP) 121C EXPOSE NEWLY ACCRETED MATTER TO 1/2 BURNING 122 X2NEW = XH2NACC/FSTEP*EXP(-RATE*TSTEP) 123 XH2NA = XH2NA + X2NEW 124 R0 = R0 + DELRATE 125 END DO 126C NOW PERFORM A WEIGHTED SUM OF THE ORIGINAL 127C AND NEWLY ACCRETED MATTER 128 DXH2O = XH2NO - XH2 129 DXH2A = XH2NA - COMPACC(12) 130 XH2N = (XH2NO*SUMM+XH2NA*FMASSACC)/(SUMM+FMASSACC) 131 XN = X + 0.5D0*(DXH2O*SUMM+DXH2A*FMASSACC)/(SUMM+FMASSACC) 132 XHE3N = XHE3 -1.5D0*(DXH2O*SUMM+DXH2A*FMASSACC)/(SUMM+FMASSACC) 133 ELSE 134C INCREMENT H,D,HE3 WITHOUT MASS ACCRETION 135 DXH2 = XH2NO - XH2 136 XH2N = XH2NO 137 XN = X + 0.5D0*DXH2 138 XHE3N = XHE3 - 1.5D0*DXH2 139 ENDIF 140C IF(IEND.EQ.M)THEN 141C WRITE(*,911)IBEGIN,IEND,DRATE(IBEGIN),RATE0,RATE1, 142C * XH2,XH2N,XHE3,XHE3N 143C WRITE(*,*)DXH2O,DXH2A,FMASSACC 144C ENDIF 145C 911 FORMAT(1X,'DBURN: MIX D BURNING',2I5,1P7E12.3) 146 DO I = IBEGIN,IEND 147 HCOMP(1,I) = XN 148 HCOMP(4,I) = XHE3N 149 HCOMP(12,I) = XH2N 150 END DO 151 9999 9999 CONTINUE 152 RETURN 153 END

1"""Python version of dburn.f (YREC 7be4133): deuterium burning over a 2step in one shell or one mixed zone.""" 3 4import math 5import numpy as np 6 7from pyrec import deuter_mod 8from pyrec import nmlpar_mod 9from pyrec import oldmod_mod 10 11 12def dburn(ibegin: int, iend: int, m: int, hs2: np.ndarray, hcomp: np.ndarray, dt: float) -> None: 13 """Burn deuterium over the step in shells ibegin..iend (one shell, or a 14 zone mixed to its mass-weighted means), with the rate rising linearly 15 from its start-of-step to its end-of-step value in substeps; update H 16 and 3He, and at the surface with accretion burn the accreted D too. 17 18 Translates: dburn.f#DBURN (7be4133) 19 20 Arguments: ibegin, iend, the shells; m, the number of shells; hs2, the 21 shell masses; hcomp[species, shell], changed in place (species 1, 4, 22 12: H, 3He, D); dt, the step (Gyr). 23 24 Returns nothing. 25 26 """ 27 28 # Module data (Fortran lines 11–14). 29 hcompp = oldmod_mod.hcompp 30 drate = deuter_mod.drate 31 drate0 = deuter_mod.drate0 32 33 # Means over the region: one shell, or mass-weighted over a zone 34 # (Fortran lines 43–74). 35 if ibegin == iend: 36 summ = hs2[ibegin] 37 x = hcompp[1, ibegin] 38 xh2 = hcompp[12, ibegin] 39 xhe3 = hcompp[4, ibegin] 40 rate0 = drate0[ibegin] 41 rate1 = drate[ibegin] 42 else: 43 summ = 0.0 44 sumr0 = 0.0 45 sumr1 = 0.0 46 x = 0.0 47 xh2 = 0.0 48 xhe3 = 0.0 49 for i in range(ibegin, iend + 1): 50 summ = summ + hs2[i] 51 sumr0 = sumr0 + hs2[i] * drate0[i] 52 sumr1 = sumr1 + hs2[i] * drate[i] 53 x = x + hcompp[1, i] * hs2[i] 54 xh2 = xh2 + hcompp[12, i] * hs2[i] 55 xhe3 = xhe3 + hcompp[4, i] * hs2[i] 56 rate0 = sumr0 / summ 57 rate1 = sumr1 / summ 58 x = x / summ 59 xh2 = xh2 / summ 60 xhe3 = xhe3 / summ 61 62 # Too little deuterium (with the accreted D at the surface): zero it 63 # and return (Fortran lines 75–85). 64 fmassacc = deuter_mod.fmassacc 65 compacc = nmlpar_mod.compacc 66 if nmlpar_mod.lmdot and iend == m and nmlpar_mod.dmdt0 > 0.0: 67 xh2test = (xh2 * summ + compacc[12] * fmassacc) / (summ + fmassacc) 68 else: 69 xh2test = xh2 70 if xh2test < 1.0e-11: 71 for i in range(ibegin, iend + 1): 72 hcomp[12, i] = 0.0 73 return # label 9999 74 75 # Burn in substeps, the rate rising linearly from rate0 to rate1 76 # (Fortran lines 86–104). 77 # Question (see the page): the comment says 0.5 in the natural log per 78 # substep; the code limits rate*tstep to 0.1. 79 # Question (see the page): the factor 2 in the exponent. 80 rate = max(rate0, rate1) 81 istep = int(rate * dt / 0.1) + 1 82 fstep = float(istep) 83 tstep = dt / fstep 84 r0 = rate0 85 delrate = (rate1 - rate0) / fstep 86 xh2no = xh2 87 for i in range(1, istep + 1): 88 rate = r0 + 0.5 * delrate 89 xh2no = xh2no * math.exp(-2.0 * rate * tstep) 90 r0 = r0 + delrate 91 92 # With accretion, burn the accreted D too and take mass-weighted means; 93 # without, the H and 3He bookkeeping alone (Fortran lines 105–139). 94 # Suspected bug (see the page): unlike the test above, this branch does 95 # not check dmdt0 > 0, so under mass loss it reads an fmassacc that MDOT 96 # never set. 97 if nmlpar_mod.lmdot and iend == m: 98 # Accreted matter: earlier accretion burns fully, new accretion is 99 # exposed to half a substep (Fortran lines 114–125). 100 r0 = rate0 101 xh2nacc = compacc[12] 102 xh2na = 0.0 103 for i in range(1, istep + 1): 104 rate = r0 + 0.5 * delrate 105 xh2na = xh2na * math.exp(-2.0 * rate * tstep) 106 x2new = xh2nacc / fstep * math.exp(-rate * tstep) 107 xh2na = xh2na + x2new 108 r0 = r0 + delrate 109 110 # Mass-weighted means of the original and the accreted matter 111 # (Fortran lines 128–132). 112 dxh2o = xh2no - xh2 113 dxh2a = xh2na - compacc[12] 114 xh2n = (xh2no * summ + xh2na * fmassacc) / (summ + fmassacc) 115 xn = x + 0.5 * (dxh2o * summ + dxh2a * fmassacc) / (summ + fmassacc) 116 xhe3n = xhe3 - 1.5 * (dxh2o * summ + dxh2a * fmassacc) / (summ + fmassacc) 117 else: 118 # Without accretion (Fortran lines 134–138). 119 dxh2 = xh2no - xh2 120 xh2n = xh2no 121 xn = x + 0.5 * dxh2 122 xhe3n = xhe3 - 1.5 * dxh2 123 124 # The new H, 3He and D over the region (Fortran lines 146–150). 125 for i in range(ibegin, iend + 1): 126 hcomp[1, i] = xn 127 hcomp[4, i] = xhe3n 128 hcomp[12, i] = xh2n

dburn.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (harmless in the runs) FMASSACC is read under mass loss.

Definitions

subroutineDBURN

Change log

1. (harmless in the runs) `FMASSACC` is read under mass loss.