dburnm.f

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

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

Variables on this page (6)

compacc, csecyr, fmassacc, hcompp, json, lmdot

Purpose

The same deuterium burning as DBURN, for the substeps of a model in the middle of a time step: it is called by MIDMOD with rates interpolated to the substep and with the accreted mass scaled to the substep's share ("unlike the routine called by MIX, the timestep here is only a portion of the total model DT"). Its header still says "DBURN."

How it fits

Inputs: as DBURN, plus DRATEM, DRATEM0 (the rates at the end and start of the substep, interpolated in MIDMOD from DRATE0, DRATE with the fraction FAC2) and FAC2 itself. DT is in seconds (MIDMOD's FAC2 = DT/DELTS), converted to Gyr with DTT = DT*1.0D-9/CSECYR (CSECYR from CONSTANTS_MOD, formerly /CONST3/). The other module data are DBURN's: HCOMPP from OLDMOD_MOD, FMASSACC from DEUTER_MOD, LMDOT and COMPACC from NMLPAR_MOD.
Who calls it: MIDMOD.
Output: HCOMP (H, $^3$He, D) over the shell or zone.

The method

As in DBURN: mass-weighted means over the zone, substeps with $X_D \leftarrow X_D e^{-2r_k\Delta t/n}$ (the factor 2 as in DBURN, whose question 2 asks what it is for, dburn.f, question 2), then $\Delta X = \tfrac12\Delta X_D$ and $\Delta X_3 = -\tfrac32\Delta X_D$. The differences from DBURN:

the threshold for zeroing deuterium is $10^{-14}$ (DBURN: $10^{-11}$) and does not include the accreted D;
with accretion, the substep's share of the accreted matter is $f_a$ = FAC2*FMASSACC, in grams, since FMASSACC is the accreted mass (MDOT sets it to DELM); with $M$ = SUMM the region's mass (for a single shell, its own mass), the new mean is weighted by mass, as in DBURN: $$X_D = \frac{X_{D,\mathrm{old}}\,M + f_a X_{D,\mathrm{acc}}}{M + f_a}.$$

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) `FMASSACC` is read under mass loss.

(line 109 ↗; lines 140–145 ↗) As in DBURN (dburn.f, suspected bug 1): The accretion branch does not test DMDT0 > 0, so in a run that only loses mass it blends with an FMASSACC that MDOT never set (mdot.f, line 320 ↗, deuter_mod.f, line 30 ↗); with the usual zero start value the result is the same as without accretion. Only with LMDOT (off by default).

A possible fix:

dburnm.f

@@ lines 106–112 @@
          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 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): 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). DBURNM runs only in MIDMOD, which only GETW calls, in rotating runs, and the two made-up mass-loss runs (YREC's 3 $M_\odot$ main-sequence and 1 $M_\odot$ giant-branch cases with LMDOT and LREIMER on) do not rotate. Reaching it needs LMDOT with LROT, which no case of YREC's sets. In DBURN, where the same branch was taken under mass loss (dburn.f, suspected bug 1), FMASSACC was 0 and the result unchanged.

Definitions

subroutineDBURNM

lines 4–170 · GitHub

Deuterium burning over a substep in MIDMOD, with interpolated rates and the substep's share of accretion; otherwise as DBURN.

Uses modules CONSTANTS_MOD, DEUTER_MOD, DIMENSIONS_MOD, NMLPAR_MOD, OLDMOD_MOD

Called by MIDMOD ×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. The accretion blend is weighted by mass, $(X M + X_\mathrm{acc} f_a)/(M + f_a)$, as in DBURN; before, it divided by $1 + f_a$ with $f_a$ in grams, so any accretion replaced the region's D, H, and $^3$He by the accreted values. For a single shell, the region mass SUMM is now set to the shell's own mass (before, a value left from an earlier call). This settles former question 2.

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

1"""Python version of dburnm.f (YREC 7be4133): deuterium burning over a 2substep of a model in the middle of a time step (MIDMOD).""" 3 4import math 5import numpy as np 6 7from pyrec import constants_mod 8from pyrec import deuter_mod 9from pyrec import nmlpar_mod 10from pyrec import oldmod_mod 11 12 13def dburnm(ibegin: int, iend: int, m: int, hs2: np.ndarray, hcomp: np.ndarray, dt: float, dratem: np.ndarray, 14 dratem0: np.ndarray, fac2: float) -> None: 15 """As DBURN, for a substep: burn deuterium in shells ibegin..iend with 16 rates interpolated to the substep, and the substep's share fac2 of the 17 accreted matter at the surface. 18 19 Translates: dburnm.f#DBURNM (7be4133) 20 21 Arguments: ibegin, iend, the shells; m, the number of shells; hs2, the 22 shell masses; hcomp[species, shell], changed in place (H, 3He, D); dt, 23 the substep (s); dratem, dratem0, the rates at the end and start of the 24 substep; fac2, the substep's fraction of the step. 25 26 Returns nothing. 27 28 """ 29 30 # The substep in Gyr (Fortran line 50). 31 hcompp = oldmod_mod.hcompp 32 dtt = dt * 1.0e-9 / constants_mod.csecyr 33 34 # Means over the region: one shell, or mass-weighted over a zone 35 # (Fortran lines 51–82). 36 if ibegin == iend: 37 summ = hs2[ibegin] 38 x = hcompp[1, ibegin] 39 xh2 = hcompp[12, ibegin] 40 xhe3 = hcompp[4, ibegin] 41 rate0 = dratem0[ibegin] 42 rate1 = dratem[ibegin] 43 else: 44 summ = 0.0 45 sumr0 = 0.0 46 sumr1 = 0.0 47 x = 0.0 48 xh2 = 0.0 49 xhe3 = 0.0 50 for i in range(ibegin, iend + 1): 51 summ = summ + hs2[i] 52 sumr0 = sumr0 + hs2[i] * dratem0[i] 53 sumr1 = sumr1 + hs2[i] * dratem[i] 54 x = x + hcompp[1, i] * hs2[i] 55 xh2 = xh2 + hcompp[12, i] * hs2[i] 56 xhe3 = xhe3 + hcompp[4, i] * hs2[i] 57 rate0 = sumr0 / summ 58 rate1 = sumr1 / summ 59 x = x / summ 60 xh2 = xh2 / summ 61 xhe3 = xhe3 / summ 62 63 # Too little deuterium: zero it and return (Fortran lines 83–88). 64 # Question (see the page): 1e-14 here, 1e-11 and the accreted D in DBURN. 65 if xh2 < 1.0e-14: 66 for i in range(ibegin, iend + 1): 67 hcomp[12, i] = 0.0 68 return # label 9999 69 70 # Burn in substeps, the rate rising linearly from rate0 to rate1 71 # (Fortran lines 89–107). 72 rate = max(rate0, rate1) 73 istep = int(rate * dtt / 0.1) + 1 74 fstep = float(istep) 75 tstep = dtt / fstep 76 r0 = rate0 77 delrate = (rate1 - rate0) / fstep 78 xh2no = xh2 79 for i in range(1, istep + 1): 80 rate = r0 + 0.5 * delrate 81 xh2no = xh2no * math.exp(-2.0 * rate * tstep) 82 r0 = r0 + delrate 83 84 # With accretion, burn the substep's share of the accreted D too and 85 # take mass-weighted means; without, the H and 3He bookkeeping alone 86 # (Fortran lines 108–162). 87 # Suspected bug (see the page): no dmdt0 > 0 test, so under mass loss 88 # this reads an fmassacc that MDOT never set. 89 compacc = nmlpar_mod.compacc 90 if nmlpar_mod.lmdot and iend == m: 91 # Accreted matter (Fortran lines 117–128). 92 r0 = rate0 93 xh2nacc = compacc[12] 94 xh2na = 0.0 95 for i in range(1, istep + 1): 96 rate = r0 + 0.5 * delrate 97 xh2na = xh2na * math.exp(-2.0 * rate * tstep) 98 x2new = xh2nacc / fstep * math.exp(-rate * tstep) 99 xh2na = xh2na + x2new 100 r0 = r0 + delrate 101 102 # Mass-weighted means with the substep's share of the accreted 103 # mass (Fortran lines 140–145). 104 fmassa = fac2 * deuter_mod.fmassacc 105 dxh2o = xh2no - xh2 106 dxh2a = xh2na - compacc[12] 107 xh2n = (xh2no * summ + xh2na * fmassa) / (summ + fmassa) 108 xn = x + 0.5 * (dxh2o * summ + dxh2a * fmassa) / (summ + fmassa) 109 xhe3n = xhe3 - 1.5 * (dxh2o * summ + dxh2a * fmassa) / (summ + fmassa) 110 else: 111 # Without accretion (Fortran lines 153–156). 112 dxh2 = xh2no - xh2 113 xh2n = xh2no 114 xn = x + 0.5 * dxh2 115 xhe3n = xhe3 - 1.5 * dxh2 116 117 # The new H, 3He and D over the region (Fortran lines 163–167). 118 for i in range(ibegin, iend + 1): 119 hcomp[1, i] = xn 120 hcomp[4, i] = xhe3n 121 hcomp[12, i] = xh2n

dburnm.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (not reached) FMASSACC is read under mass loss.

Definitions

subroutineDBURNM

Change log

1. (not reached) `FMASSACC` is read under mass loss.