microdiff_run.f

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

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

Variables on this page (6)

c4pi, grtol, iowr, ishort, json, niter_gs

Purpose

Diffuses one species over one time step on the equally spaced grid, for MICRODIFF, and returns the change in its abundance. The comments describe the diffusion equation as two terms: one with $d\ln P/dr$ and the abundance gradient (solved explicitly, by the two-step Lax–Wendroff scheme of Numerical Recipes), and one with the second derivative of the abundance (solved implicitly).

How it fits

Inputs: the equal grid and midpoint variables from MICRODIFF_MTE: mass EM, radius ER, density, temperature, $d\ln P/dr$, $\nabla$ and the three-species abundances ESPEC(3,·); the spacing DR; DT and STOT in Bahcall–Loeb units; and the species' weight THEA, charge THEZ and column J (1 for hydrogen, 3 for metals or a light element).
Who calls it: only MICRODIFF, once per diffused species.
What it calls: MICRODIFF_COD (the coefficients, twice), LAX_WENDROF1 and LAX_WENDROF2 (the explicit term), GET_IMP_DIFFCO (the matrix) and TRIDIAG_GS (its solution).
Output: ESPEC(J,·) is overwritten with the change $X^\mathrm{new}_i - X^\mathrm{old}_i$; EDXDR (the hydrogen concentration gradient, filled on the hydrogen call) is kept for later calls. A line per iteration goes to ISHORT, and a warning to IOWR and ISHORT if the iteration does not converge.
Modules: NMLPAR_MOD (GRTOL, NITER_GS), CONSTANTS_MOD ($4\pi$, C4PI), IOUNITS_MOD (ISHORT, IOWR) and DIMENSIONS_MOD (JSON). The unused /THEAGE/ of rp_1 is no longer declared.

The method

Let $X_i$ be the abundance at grid point $i$ (radius $r_i$, mass $m_i$), and $X_{i+1/2}$ at the midpoints.

Coefficients at the points. MICRODIFF_COD gives the advective coefficient $C_{1,i}$ (ECOD1) and the diffusion coefficient $C_{2,i}$ (ECOD2).
Lax–Wendroff, first half step (LAX_WENDROF1), the change at the midpoints: $$\Delta X_{i+1/2} = \frac{4\pi\,\Delta t}{2}\;\frac{C_{1,i+1} - C_{1,i}}{m_{i+1} - m_i},$$ where the first cell takes $C_{1,1} = 0$ and $m_1 = 0$, and the last cell $C_{1,N} = 0$ and $m_N =$ STOT, so the first and last cells reach to the centre and the surface.
Coefficients at the midpoints, from the half-step abundances $X_{i+1/2} + \Delta X_{i+1/2}$.
Lax–Wendroff, full step (LAX_WENDROF2), at the points: $$X'_i = X_i + 4\pi\,\Delta t\;\frac{C_{1,i+1/2} - C_{1,i-1/2}}{m_{i+1/2} - m_{i-1/2}},$$ with $C_{1,1/2} = m_{1/2} = 0$ at the first point, and $C_{1,N+1/2} = 0$, $m_{N+1/2} =$ STOT at the last.
The implicit term. With $\alpha_i = 4\pi\,\Delta t / (\Delta r\, \Delta m_i)$ (ALPHA, where $\Delta m_i = m_{i+1/2} - m_{i-1/2}$, from the centre for $i=1$ and to STOT for $i=N$), GET_IMP_DIFFCO builds and TRIDIAG_GS solves $$X_i - \alpha_i\left[C_{2,i+1/2}\,(X_{i+1}-X_i) - C_{2,i-1/2}\,(X_i - X_{i-1})\right] = X'_i .$$ The loop repeats until the largest change between iterations is below GRTOL, at most NITER_GS times.

Question 1 (lines 160–174 ↗; get_imp_diffco.f, line 20 ↗). An iteration that cannot change anything. Inside the loop, GET_IMP_DIFFCO corrects ECOD2_H by EQCOD2X_H * ED_H, but EQCOD2X_H is set to zero just before ("only used in GRSETT … retained for legacy"). So the matrix and the right-hand side ED_PRIME are the same in every iteration, and the second iteration always reproduces the first (its change is zero, below GRTOL). Is the loop (and the convergence test on GRTOL) still meaningful, or should the coefficients be recomputed from the new abundances?

The comments name the elements of an older use. They speak of "the original run of light element mass fractions," "the new run of LI," and "the corrections to the hydrogen abundance" (line 82 ↗; line 172 ↗; lines 176–177 ↗), but ED is loaded from ESPEC(J,·) (line 80 ↗) for whichever species J the caller passes, so they could say "the diffused species."
Result: $\Delta X_i = X^\mathrm{new}_i - X_i$ in ESPEC(J,·).

Design

The change rather than the new abundance is returned, as the comments explain, so that the interpolation back to the model grid adds less error. The first and last cells hold the masses of the convective core and envelope, whose changes MICRODIFF_ETM then applies to every point inside them.

From YREC's documentation

Physics Options: GRTOL (default 1.0D-8) is the numerical tolerance on the abundance changes from diffusion, and NITER_GS (10) the largest number of iterations of the solver.

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) With `NITER_GS = 0`, the implicit term is skipped and the warning prints old values.

(line 152 ↗; lines 178–196 ↗) With NITER_GS = 0 the loop never runs: The implicit (diffusive) term is not applied at all, and the warning "MICRODIFF FAILED TO CONVERGE" prints DXMAX and IMAX left from an earlier call (all locals are SAVEd). Only with NITER_GS = 0 (the default is 10).

A possible fix: reject NITER_GS < 1 where the namelist is read, or set DXMAX = 0.0D0 and IMAX = 0 before the loop so that the message is at least true.

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 loop was entered at every settling call in the seven runs with settling (testsuite/Test_solarGS98_dif, testsuite/Test_solarAAG21_dif, testsuite/Test_solarM22M_difkur, testsuite/Test_solarGS98_difrotmixper18p4d1p0, testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, testsuite/Test_m1p4feh+0p0GS98_dif_TAMS, and testsuite/Test_solarAAG21_dif with a custom mixture), up to 18,515 times per run, always with NITER_GS = 10. Not reached: It needs NITER_GS = 0, and every case sets 10. No fix was run, since the page's fixes only stop the run or correct the message. Intent: none found.

2. (harmless in the runs) The convergence test is one-sided.

(lines 178–190 ↗) DXMAX takes the largest signed change ED(I)-ED_P(I), starting from zero, so negative changes are never seen, in the test or in the printed line. Correct today, since the second iteration reproduces the first (question 1); it would end the loop too early if the iterations changed anything.

A possible fix:

microdiff_run.f

@@ lines 177–184 @@
 C  ENOUGH TO EXIT.
          DXMAX = 0.0D0
          IMAX = 0
+C 11/09/26 KC via Claude
+C code changed to test the size of the change, not its signed value.
          DO I = 1,NPT
-            DX = ED(I)-ED_P(I)
+            DX = ABS(ED(I)-ED_P(I))
             IF(DX.GT.DXMAX)THEN
                DXMAX = DX
                IMAX = 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 test ran in all seven runs with settling, 4,496 to 22,220 times per run, never more than two iterations a call. In 14 to 259 first iterations per run, the largest signed change was below GRTOL while the largest change in size, $2.0\times10^{-7}$ to $1.87\times10^{-6}$, was not, so the loop ended one iteration early. With the fix, testsuite/Test_m1p4feh+0p0GS98_dif_TAMS and testsuite/Test_solarM22M_difkur give identical tracks and stored models; the .short file differs beyond 2% (presumably its iteration lines; not checked). This agrees with the above: The second iteration repeats the first. Intent: none found.

3. (cannot be reached) The midpoint change is shifted by one point.

(lines 155–157 ↗) ED_H(I) = 0.5*(ED(I)+ED(I-1)-…) for I = 2..NPT assigns to midpoint I the mean of points I-1 and I, while in MICRODIFF_MTE midpoint I lies between points I and I+1 (microdiff_mte.f, lines 85–93 ↗); midpoint 1 keeps its value from LAX_WENDROF1. No effect today, since EQCOD2X_H is zero (question 1), but it would matter if EQCOD2X_H were used.

A possible fix:

microdiff_run.f

@@ lines 152–159 @@
       DO 100 ITER=1,NITER_GS
 C  FIND CHANGE IN D AT THE ZONE MIDPOINTS, GIVEN CHANGE IN D AT
 C  THE ZONE CENTERS.
-         DO I = 2,NPT
-            ED_H(I) = 0.5D0*(ED(I)+ED(I-1)-ED_P(I)-ED_P(I-1))
+C 11/09/26 KC via Claude
+C code changed to put the change between points I and I+1 at
+C midpoint I, as MICRODIFF_MTE places the midpoints.
+         DO I = 1,NPT-1
+            ED_H(I) = 0.5D0*(ED(I+1)+ED(I)-ED_P(I+1)-ED_P(I))
          ENDDO
 C  STORE CURRENT RUN OF ABUNDANCES IN VECTOR ED_P. THE ITERATION
 C  LOOP IS COMPLETED ONCE ED - ED_P IS EVERYWHERE LESS THAN THE TOLERANCE

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 midpoint changes were computed 4,496 to 22,220 times per run in all seven runs with settling, and EQCOD2X_H was zero every time. It cannot matter at 7be4133, since EQCOD2X_H is set to zero just before every use (lines 162–165 ↗); what remains is fragile code. The fix was not run. Intent: The comment there says EQCOD2X_H "IS ONLY USED IN GRSETT (OLD SETTLING ROUTINE), BUT RETAINED FOR LEGACY" (lines 160–161 ↗).

Definitions

subroutineMICRODIFF_RUN

lines 5–204 · GitHub

One species, one step: explicit two-step Lax–Wendroff for the advective term, then an implicit tridiagonal solve for the diffusive term (equations above). Returns the change in abundance.

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD

Called by MICRODIFF ×3

Flowchart of its calls, in order (6)

MICRODIFF_COD 101→LAX_WENDROF1 110→MICRODIFF_COD 124→LAX_WENDROF2 133

for ITER=1,NITER_GS (DO 100), lines 152–191

GET_IMP_DIFFCO 170→TRIDIAG_GS 174

Change log

The changes since rp_1

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

1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2C MICRODIFF_RUN 3CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4C SUBROUTINE MICRODIFF_RUN(NLIGHT,CSLIGHT,DR,DT,STOT,NPT,EM,ER,EDEN,ETEM, ! KC 2025-05-31 5 SUBROUTINE MICRODIFF_RUN(DR,DT,STOT,NPT,EM,ER,EDEN,ETEM, 6 * * EHQPR,EDELR,ESPEC,EDXDR,EM_H,ER_H,EDEN_H,ETEM_H,EHQPR_H, 7 * * EDELR_H,ESPEC_H,EDXDR_H,THEA,THEZ,J) 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/NMLPAR/ into NMLPAR_MOD. 14 USE NMLPAR_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 17 USE CONSTANTS_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 20 USE IOUNITS_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 A,ALPHA,B,C,DR,DT,DX,DXMAX,ECOD1,ECOD1_H,ECOD2,ECOD2_H,ED, 26 * * EDELR,EDELR_H,EDEN,EDEN_H,EDXDR,EDXDR_H,ED_H,ED_ORIG, 27 * * ED_ORIG_H,ED_P,ED_PRIME,EHQPR,EHQPR_H,EM,EM_H,EQCOD2X_H,ER, 28 * * ER_H,ESPEC,ESPEC_H,ETEM,ETEM_H,FAC,STOT,THEA,THEZ 29 INTEGER I,IMAX,ITER,J,NPT,NPT2 30 LOGICAL*4 LDOLI 31 32C 10/12/26 MHP via Claude 33C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 34C namelist variables that the physics routines read from 27 COMMON 35C blocks which also held derived or working members. Those blocks keep 36C only their non-namelist members; five of them -- NEWMX, VMULT, 37C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 38C Derivation-only namelist variables from the same blocks went to 39C COMMON/NMLBYPASS/ instead, not here. 40C 41C Declared identically in all 121 files that hold any of these members. 42C Generated from one member list, so every copy is byte-identical; the 43C linker's 'size of symbol changed' warning is the check that they are. 44C 10/13/26 MHP via Claude 45C code changed to drop COMMON/THEAGE/, which this routine declared but 46C never used. DAGE, the stellar age, is passed explicitly elsewhere. 47C MHP 3/94 ADDED METAL DIFFUSION 48C 9/25/26 MHP via Claude 49C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 50C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 51C 9/25/26 MHP via Claude 52C code changed to consolidate all I/O logical unit numbers, previously 53C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 54C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 55C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 56C COMMON/IOUNITS/ 57C DIMENSION HCOMP(15,JSON),HQPR(JSON),HR(JSON),HS1(JSON),HT(JSON), 58C * HD(JSON),HRU(JSON),HTU(JSON),LC(JSON) 59 DIMENSION EM(JSON),ER(JSON),EDEN(JSON),ETEM(JSON),EHQPR(JSON), 60 * * EDELR(JSON),ESPEC(3,JSON),EDXDR(JSON),EM_H(JSON),ER_H(JSON), 61 * * EDEN_H(JSON),ETEM_H(JSON),EHQPR_H(JSON),EDELR_H(JSON), 62 * * ESPEC_H(3,JSON),EDXDR_H(JSON) 63 DIMENSION ED(JSON),ED_H(JSON),ED_ORIG(JSON),ED_ORIG_H(JSON), 64 * * ED_P(JSON),ED_PRIME(JSON),ECOD1(JSON), ! ED_HALF(JSON) 65 * * ECOD1_H(JSON),ECOD2(JSON),ECOD2_H(JSON), 66 * * EQCOD2X_H(JSON) 67 DIMENSION ALPHA(JSON),A(JSON),B(JSON),C(JSON) 68 SAVE 69C 70C THIS PROGRAM EXECUTES THE DIFFUSING OF ELEMENTS. ALL OF THE RELEVANT 71C PARAMETERS ARE DEFINED AT BOTH THE GRID POINTS AND MIDPOINTS. 72C 73C CREATE DUMMY VECTORS FOR THE DIFFUSING ELEMENT. THIS ALLOWS A SINGLE 74C ELEMENT TO BE EASILY PASSED INTO SUB-ROUTINES. 75C 76 DO I=1,NPT-1 77 ED(I) = ESPEC(J,I) 78 ED_H(I) = ESPEC_H(J,I) 79 ENDDO 80 ED(NPT) = ESPEC(J,NPT) 81C 82C STORE ORIGINAL RUN OF LIGHT ELEMENT MASS FRACTIONS IN VECTOR ED_ORIG. 83C THE CHANGE IN ELEMENTS (ED - ED_ORIG) IS INTERPOLATED BACK TO THE 84C ORIGINAL GRID AT THE END OF THE ROUTINE; THE CHANGE IN ABUNDANCE 85C (RATHER THAN THE NEW RUN OF ABUNDANCE) IS USED TO MINIMIZE ERRORS 86C ARISING FROM THE INTERPOLATION. 87 DO I = 1,NPT-1 88 ED_ORIG(I) = ED(I) 89 ED_ORIG_H(I) = ED_H(I) 90 ENDDO 91 ED_ORIG(NPT) = ED(NPT) 92c 93C SET THE FLAG LDOLI = TRUE, SO THAT THE LAX-WENDROF ROUTINES KNOW TO NOT 94C USE THE OLD METAL DIFFUSION VECTORS. 95 LDOLI = .TRUE. 96c 97C NOW USE THE EQUALLY SPACED GRID TO CALUCLATE DIFFUSION COEFFICIENTS IN THE 98C METHOD OF THE THOUL ROUTINE. MUCH OF THE FOLLOWING CODE WAS TAKEN FROM 99c SETUP_LISETT.F. 100C 101 CALL MICRODIFF_COD(NPT,ESPEC,ER,EDEN,ETEM,EHQPR,EDELR,ECOD1,ECOD2, 102 * * EDXDR,THEA,THEZ,J) 103C 104C FIRST STEP OF TWO STEP LAX-WENDROFF METHOD. COMPUTE NEW ABUNDANCES AT ZONE 105C MIDPOINTS USING THE LAX SCHEME : 106C X(N+1/2,J+1/2)=1/2(X(N,J+1/2)-1/2(DT/DR)(COD1(N,J+1)-COD1(N,J)) 107C WHERE N IS THE TIME VARIABLE, J IS THE SPATIAL ONE, AND COD1 IS THE 108C DIFFUSION COEFFICIENT. 109C 110 CALL LAX_WENDROF1(DT,ECOD1,EM,NPT,STOT,ED_H,LDOLI) 111C 112C UPDATE ESPEC WITH THE NEW RUN OF ED_H, FOR THE PURPOSE OF 113C CALUCLATING NEW DIFFUSION COEFFICIENTS. 114C 115 DO I=1,NPT-1 116 ESPEC_H(J,I) = ED_ORIG_H(I) + ED_H(I) 117 ENDDO 118C 119C GET NEW DIFFUSION COEFFICIENTS AT THE ZONE MIDPOINTS, USING THE 120C PROVISIONAL SOLUTION. NOTE THAT NPT-1 MUST BE PASSED IN, AS THE 121c MIDPOINTS HAVE ONE LESS GRID ELEMENT. 122C 123 NPT2 = NPT-1 124 CALL MICRODIFF_COD(NPT2,ESPEC_H,ER_H,EDEN_H,ETEM_H,EHQPR_H, 125 * * EDELR_H,ECOD1_H,ECOD2_H,EDXDR_H,THEA,THEZ,J) 126C 127C USING THE NEW COEFFICIENTS, SOLVE FOR THE NEW RUN OF ABUNDANCES. 128C NOTE : THE SR ACTUALLY RETURNS THE *CHANGE* IN THE ABUNDANCE AS A FUNCTION 129C OF RADIUS, WHICH IS APPLIED TO THE ORIGINAL RUN OF ELEMENTS, RATHER THAN 130C COMPUTING A NEW RUN OF ABUNDANCE AND TRANSFORMING IT BACK. THIS IS DONE 131C TO MINIMIZE ERRORS FROM THE INTERPOLATION. 132C 133 CALL LAX_WENDROF2(DT,ECOD1_H,EM_H,ED,NPT,STOT,LDOLI) 134C 135C NOW IMPLICITLY SOLVE FOR THE SECOND TERM (INVOLVING THE SECOND 136C DERIVATIVE OF THE COMPOSITION GRADIENT). FIRST DEFINE VECTOR ED_P, 137C WHICH IS THE ABUNDANCE AT THE START OF THE TIMESTEP, AND ELI_PRIME, 138C THE ABUNDANCE ONE WOULD HAVE WITHOUT THE SECOND TERM IN THE DIFFUSION 139C EQUATION. 140 DO I = 1,NPT 141 ED_P(I) = ED_ORIG(I) 142 ED_PRIME(I) = ED(I) 143 ENDDO 144C ALPHA IS THE NUMERICAL FACTOR IN FRONT OF THE DIFFUSION COEFFICIENTS. 145 FAC = C4PI*DT/DR 146 ALPHA(1) = FAC/EM_H(1) 147 DO I = 2,NPT-1 148 ALPHA(I) = FAC/(EM_H(I)-EM_H(I-1)) 149 ENDDO 150 ALPHA(NPT) = FAC/(STOT-EM_H(NPT-1)) 151C START ITERATION LOOP FOR THE NEW RUN OF ABUNDANCES. 152 DO 100 ITER=1,NITER_GS 153C FIND CHANGE IN D AT THE ZONE MIDPOINTS, GIVEN CHANGE IN D AT 154C THE ZONE CENTERS. 155 DO I = 2,NPT 156 ED_H(I) = 0.5D0*(ED(I)+ED(I-1)-ED_P(I)-ED_P(I-1)) 157 ENDDO 158C STORE CURRENT RUN OF ABUNDANCES IN VECTOR ED_P. THE ITERATION 159C LOOP IS COMPLETED ONCE ED - ED_P IS EVERYWHERE LESS THAN THE TOLERANCE 160C GRTOL. ALSO ZERO OUT EQCOD2X_H, WHICH IS ONLY USED IN GRSETT (OLD SETTLING 161C ROUTINE), BUT RETAINED FOR LEGACY. 162 DO I = 1,NPT-1 163 ED_P(I) = ED(I) 164 EQCOD2X_H(I) = 0.0 165 ENDDO 166 ED_P(NPT) = ED(NPT) 167C GET NEW DIFFUSION COEFFICIENTS, TAKING INTO ACCOUNT THE CHANGE IN D 168C FROM THE PREVIOUS ITERATION. 169C 170 CALL GET_IMP_DIFFCO(ALPHA,ECOD2_H,ED_H,EQCOD2X_H,A,B,C,NPT) 171C 172C SOLVE THE TRIDIAGONAL MATRIX SYSTEM FOR THE NEW RUN OF LI. 173C 174 CALL TRIDIAG_GS(A,B,C,ED_PRIME,NPT,ED) 175C 176C CHECK TO SEE IF THE CORRECTIONS TO THE HYDROGEN ABUNDANCE ARE SMALL 177C ENOUGH TO EXIT. 178 DXMAX = 0.0D0 179 IMAX = 0 180 DO I = 1,NPT 181 DX = ED(I)-ED_P(I) 182 IF(DX.GT.DXMAX)THEN 183 DXMAX = DX 184 IMAX = I 185 ENDIF 186 ENDDO 187 WRITE(ISHORT,90)ITER,DXMAX,IMAX 188 90 90 FORMAT(1X,'ITERATION ',I3,' DXMAX ',1PE10.2,' IMAX ',I4) 189C EXIT ITERATION LOOP IF SYSTEM HAS CONVERGED. 190 IF(DXMAX.LT.GRTOL)GOTO 120 191 100 100 CONTINUE 192 WRITE(IOWR,110)GRTOL,NITER_GS,DXMAX,IMAX 193 WRITE(ISHORT,110)GRTOL,NITER_GS,DXMAX,IMAX 194 110 110 FORMAT(1X,'MICRODIFF FAILED TO CONVERGE TO WITHIN ',1PE9.3,' IN ',I3, 195 * * 'ITERATIONS'/1X,'LAST ITERATION CHANGE IN D ',1PE9.3, 196 * * ' IN EQUALLY SPACED SHELL ',I5) 197 120 120 CONTINUE 198C 199C STORE THE RUN OF CHANGES IN THE DIFFUSED ELEMENT. 200 DO I = 1,NPT 201 ESPEC(J,I) = ED(I)-ED_ORIG(I) 202 ENDDO 203 RETURN 204 END

1"""Python version of microdiff_run.f (YREC 7be4133): one species diffused 2over one step on the equally spaced grid, for MICRODIFF.""" 3 4import numpy as np 5 6from pyrec import constants_mod 7from pyrec import dimensions_mod 8from pyrec import iounits_mod 9from pyrec import nmlpar_mod 10from pyrec.get_imp_diffco import get_imp_diffco 11from pyrec.lax_wendrof1 import lax_wendrof1 12from pyrec.lax_wendrof2 import lax_wendrof2 13from pyrec.microdiff_cod import microdiff_cod 14from pyrec.tridiag_gs import tridiag_gs 15 16JSON = dimensions_mod.json 17 18# The Fortran's local arrays (SAVE). 19ed = np.zeros(JSON + 1) 20ed_h = np.zeros(JSON + 1) 21ed_orig = np.zeros(JSON + 1) 22ed_orig_h = np.zeros(JSON + 1) 23ed_p = np.zeros(JSON + 1) 24ed_prime = np.zeros(JSON + 1) 25ecod1 = np.zeros(JSON + 1) 26ecod1_h = np.zeros(JSON + 1) 27ecod2 = np.zeros(JSON + 1) 28ecod2_h = np.zeros(JSON + 1) 29eqcod2x_h = np.zeros(JSON + 1) 30alpha = np.zeros(JSON + 1) 31a = np.zeros(JSON + 1) 32b = np.zeros(JSON + 1) 33c = np.zeros(JSON + 1) 34 35 36def microdiff_run(dr: float, dt: float, stot: float, npt: int, em: np.ndarray, er: np.ndarray, eden: np.ndarray, 37 etem: np.ndarray, ehqpr: np.ndarray, edelr: np.ndarray, espec: np.ndarray, edxdr: np.ndarray, 38 em_h: np.ndarray, er_h: np.ndarray, eden_h: np.ndarray, etem_h: np.ndarray, ehqpr_h: np.ndarray, 39 edelr_h: np.ndarray, espec_h: np.ndarray, edxdr_h: np.ndarray, thea: float, thez: float, 40 j: int) -> None: 41 """Diffuse species column j over one step on the equal grid: the 42 advective term by the two-step Lax-Wendroff scheme, then the diffusive 43 term implicitly (a tridiagonal solve, iterated up to NITER_GS times to 44 the tolerance GRTOL). Replaces espec[j, 1..npt] by the change in the 45 abundance. 46 47 Translates: microdiff_run.f#MICRODIFF_RUN (7be4133) 48 49 Arguments: the spacing dr; dt and stot in Bahcall and Loeb's units; 50 npt grid points; the grid variables em, er, eden, etem, ehqpr, edelr 51 and their midpoint versions (_h, npt - 1 of them); the abundances 52 espec[k, i] and espec_h[k, i] (k = 1..3); the hydrogen gradient edxdr, 53 edxdr_h (set by MICRODIFF_COD on the hydrogen call); the species' 54 weight thea, charge thez and column j. Arrays have index 0 unused. 55 56 Sets in place: espec[j, i] (the change), espec_h[j, i] (the half-step 57 abundances), edxdr, edxdr_h. Returns nothing. 58 59 """ 60 61 # Dummy vectors for the diffusing element (Fortran lines 73–80). 62 for i in range(1, npt): 63 ed[i] = espec[j, i] 64 ed_h[i] = espec_h[j, i] 65 ed[npt] = espec[j, npt] 66 67 # Its original run (Fortran lines 82–91). 68 # Question (see the page): the comments name light elements. 69 for i in range(1, npt): 70 ed_orig[i] = ed[i] 71 ed_orig_h[i] = ed_h[i] 72 ed_orig[npt] = ed[npt] 73 74 # LDOLI: the Lax-Wendroff routines do not use the old metal vectors 75 # (Fortran lines 93–95). 76 ldoli = True 77 78 # Coefficients at the grid points (Fortran lines 97–102). 79 microdiff_cod(npt, espec, er, eden, etem, ehqpr, edelr, ecod1, ecod2, 80 edxdr, thea, thez, j) 81 82 # Lax-Wendroff, first half step: the change at the midpoints, and the 83 # half-step abundances (Fortran lines 104–117). 84 lax_wendrof1(dt, ecod1, em, npt, stot, ed_h, ldoli) 85 for i in range(1, npt): 86 espec_h[j, i] = ed_orig_h[i] + ed_h[i] 87 88 # Coefficients at the npt - 1 midpoints, from the provisional solution 89 # (Fortran lines 119–125). 90 npt2 = npt - 1 91 microdiff_cod(npt2, espec_h, er_h, eden_h, etem_h, ehqpr_h, edelr_h, 92 ecod1_h, ecod2_h, edxdr_h, thea, thez, j) 93 94 # Lax-Wendroff, full step at the points (Fortran lines 127–133). 95 lax_wendrof2(dt, ecod1_h, em_h, ed, npt, stot, ldoli) 96 97 # The implicit term: ed_p is the abundance at the start of the step, 98 # ed_prime the abundance without the second term (Fortran lines 99 # 135–143). 100 for i in range(1, npt + 1): 101 ed_p[i] = ed_orig[i] 102 ed_prime[i] = ed[i] 103 104 # ALPHA, the factor in front of the diffusion coefficients (Fortran 105 # lines 144–150). 106 fac = constants_mod.c4pi * dt / dr 107 alpha[1] = fac / em_h[1] 108 for i in range(2, npt): 109 alpha[i] = fac / (em_h[i] - em_h[i - 1]) 110 alpha[npt] = fac / (stot - em_h[npt - 1]) 111 112 # The iteration for the new run of abundances (Fortran lines 151–197). 113 # Question (see the page): every iteration solves the same system. 114 # Suspected bug (see the page): with NITER_GS = 0 the Fortran prints 115 # DXMAX and IMAX left from an earlier call; here they are 0. 116 ishort = iounits_mod.ishort 117 converged = False 118 dxmax = 0.0 119 imax = 0 120 for it in range(1, nmlpar_mod.niter_gs + 1): # loop 100 121 # The change at the midpoints (Fortran lines 153–157). 122 # Suspected bug (see the page): midpoint i is taken between points 123 # i-1 and i here, between i and i+1 in MICRODIFF_MTE. 124 for i in range(2, npt + 1): 125 ed_h[i] = 0.5 * (ed[i] + ed[i - 1] - ed_p[i] - ed_p[i - 1]) 126 127 # Store the current run in ed_p, and zero eqcod2x_h, used only by 128 # GRSETT (Fortran lines 158–166). 129 for i in range(1, npt): 130 ed_p[i] = ed[i] 131 eqcod2x_h[i] = 0.0 132 ed_p[npt] = ed[npt] 133 134 # The tridiagonal system and its solution (Fortran lines 167–174). 135 get_imp_diffco(alpha, ecod2_h, ed_h, eqcod2x_h, a, b, c, npt) 136 tridiag_gs(a, b, c, ed_prime, npt, ed) 137 138 # The convergence test (Fortran lines 176–190). 139 # Suspected bug (see the page): the largest signed change, not 140 # |change|. 141 dxmax = 0.0 142 imax = 0 143 for i in range(1, npt + 1): 144 dx = ed[i] - ed_p[i] 145 if dx > dxmax: 146 dxmax = dx 147 imax = i 148 iounits_mod.write(ishort, f" ITERATION {it:3d} DXMAX " 149 f"{dxmax:10.2E} IMAX {imax:4d}") # format 90 150 if dxmax < nmlpar_mod.grtol: 151 converged = True 152 break # GOTO 120 153 154 # Not converged (Fortran lines 192–196). 155 if not converged: 156 text = (f" MICRODIFF FAILED TO CONVERGE TO WITHIN " 157 f"{nmlpar_mod.grtol:9.3E} IN {nmlpar_mod.niter_gs:3d}" 158 f"ITERATIONS\n LAST ITERATION CHANGE IN D {dxmax:9.3E}" 159 f" IN EQUALLY SPACED SHELL {imax:5d}") # format 110 160 iounits_mod.write(iounits_mod.iowr, text) 161 iounits_mod.write(ishort, text) 162 163 # Label 120: the change in the diffused element (Fortran lines 164 # 199–202). 165 for i in range(1, npt + 1): 166 espec[j, i] = ed[i] - ed_orig[i]

microdiff_run.f

Purpose

How it fits

The method

Design

From YREC's documentation

Suspected bugs, checked by running the code

1. (never occurred) With NITER_GS = 0, the implicit term is skipped and the warning prints old values.

2. (harmless in the runs) The convergence test is one-sided.

3. (cannot be reached) The midpoint change is shifted by one point.

Definitions

subroutineMICRODIFF_RUN

Change log

1. (never occurred) With `NITER_GS = 0`, the implicit term is skipped and the warning prints old values.