grsett.f

GRSETT computes the gravitational settling of helium (and, with LDIFZ, of the metals) over one time step, "using the formalism of Bahcall and Loeb 1989" (its comment). The code solves for the divergence of a flux, in mass (as LAX_WENDROF1 differences it), $$\frac{\partial X}{\partial t} = 4\pi\,\frac{\partial}{\partial m}\Bigl[D_1 + D_2\,\frac{\partial X}{\partial r}\Bigr],\qquad D_1 = f\,\frac{r^2T^{5/2}}{\ln\Lambda}\,X\,(A_P + A_T\nabla)\,\frac{d\ln P}{dr},\qquad D_2 = f\,\frac{r^2T^{5/2}}{\ln\Lambda}\,A_C,$$ where $r$ is in $R_\odot$ and $T$ in $10^7$ K (Bahcall and Loeb's units), $A_P$, $A_T$, $A_C$ are the coefficients of the pressure, temperature and concentration gradients (AP, AT/$\nabla$, AC, from Thoul et al. with LTHOUL; without it, Bahcall and Loeb's fits in $X$), $f$ the user's factor (below), and $D_1$, $D_2$ are COD1, COD2 from SETUP_GRSETT. The first term is advective (settling at a velocity proportional to $d\ln P/dr$), the second diffusive; the code's comment, a first term "that depends on DLNP/DR*DX/DR" and a second "on D^2X/DR^2," describes them only loosely. The factor $f$ is FGRY for helium and FGRZ for the metals (namelist settings, default 1), and with LTHOUL also CSTMIXING (and CSTDIFFMIX on the metals' settling term), factors added in 10/09 to "mimic mixing to reduce settling" (their $X$-derivatives are not scaled: setup_grsett.f, suspected bug 1). GRSETT treats the first term explicitly with a two-step Lax–Wendroff scheme and the second implicitly, iterating on the coefficients.

• Who calls it: MIX and ROTMIX, the mixing drivers of the evolution step.
• Inputs: DT0 (time step, s), HCOMP (mass fractions; HCOMP(1) = $X$, HCOMP(2) = $Y$, HCOMP(3) = $Z$), HQPR0 = $d\ln P/dr$, HR = $\log r$, HD = $\log\rho$ (passed on to SETUP_GRSETT, which needs it for the Coulomb logarithm), HS10 = $m$ (g, not logged), HT = $\log T$, LC (convective or not), M, STOT0. GRSETT copies DT0, STOT0, HS10 and HQPR0 into its own DT, STOT, HS1, and HQPR, which SETUP_GRSETT converts to Bahcall and Loeb's units in place and EQUAL_TO_MODEL converts back, so the caller's values are not changed.
• Output: the new HCOMP(1) and HCOMP(2), and with LDIFZ HCOMP(3) ($Z$, the CNO rows scaled with it) (via EQUAL_TO_MODEL, which uses $X+Y+Z=1$).
• Calls, in order: SETUP_GRSETT (the coefficients COD1, COD2 and their $X$-derivatives, the convection-zone edges, and whether the model is fully convective), MODEL_TO_EQUAL (to an equally spaced grid in radius), LAX_WENDROF1 and LAX_WENDROF2, then in each iteration GET_IMP_DIFFCO and TRIDIAG_GS, and at the end EQUAL_TO_MODEL.
• Modules (COMMON blocks until rp_1): NMLPAR_MOD for the settings GRTOL, NITER_GS (formerly /GRAVST/, set by PARMIN) and LDIFZ (formerly /GRAVS3/); GRAVEZ_MOD (the metals on the equal grid, with the Lax–Wendroff pair and the regridding routines); CONSTANTS_MOD; IOUNITS_MOD (the output units). The metals' coefficients COD1Z, COD2Z, QCOD1Z, QCOD2Z (formerly /GRAVSZ/) are now GRSETT's own arrays, filled by SETUP_GRSETT and passed to MODEL_TO_EQUAL.

1. Set up the coefficients; return at once for a fully convective model. Settling is computed only between the edge of a convective core (or the centre, or the first shell with $X >$ XMIN) and the base of the surface convection zone (or, if helium is too scarce there, the first shell below it with $Y >$ YMIN): IBEGIN and IEND. What lies beyond them is not on the grid; each side acts as one reservoir at the grid's end (the boundary cells of LAX_WENDROF1 hold all the mass beyond the grid), and every shell below IBEGIN (a convective core, or a hydrogen-exhausted region, which may be radiative) and above IEND gets the change found at the grid's end point (step 5); for a convective zone this keeps it homogeneous.
2. Move to $N$ (NPT) equally spaced radii $r_i$ (step $\Delta r$), with values at the points (E…) and at the midpoints (E…_H). Keep the starting $X_i^0$ (EX_ORIG).

3. Explicit term. LAX_WENDROF1 gives the half-step changes $\Delta X_{i+1/2}$; the midpoint coefficients are linearized in them, $$F_{i+1/2} \leftarrow F_{i+1/2} + \Delta X_{i+1/2}\,\frac{\partial F}{\partial X}\Big|_{i+1/2},$$ and LAX_WENDROF2 applies the full step, giving $X'_i$ (EX_PRIME).

4. Implicit term. With $\alpha_i = 4\pi\Delta t/(\Delta r\,\Delta m_i)$ (ALPHA; $\Delta m_i$ the mass between midpoints, $m_{1/2}$ at the centre and $M - m_{N-1/2}$ at the surface), iterate up to NITER_GS times: take the change at the midpoints $\tfrac12(X_i + X_{i-1} - X^p_i - X^p_{i-1})$ since the previous iterate $X^p$; let GET_IMP_DIFFCO build the tridiagonal coefficients $a_i, b_i, c_i$ from ECOD2_H corrected by that change; solve $$a_i X_{i-1} + b_i X_i + c_i X_{i+1} = X'_i$$ with TRIDIAG_GS. With $D_{i+1/2}$ = ECOD2_H(I) (the coefficient between points $i$ and $i+1$), the rows are those of an implicit (backward Euler) step: $$\bigl(1+\alpha_i(D_{i-1/2}+D_{i+1/2})\bigr)X_i - \alpha_i D_{i-1/2}X_{i-1} - \alpha_i D_{i+1/2}X_{i+1} = X'_i,$$ with $D_{1/2} = D_{N+1/2} = 0$ (no flux through the centre or the surface). Stop when $\max_i (X_i - X^p_i) <$ GRTOL, the largest signed change (suspected bug 2); print each iteration's largest change. The change at the midpoints is stored one index off from the coefficients it corrects (suspected bug 1).

Each iteration writes a line to the short output (ISHORT). The metals get the explicit term only; the implicit iteration is for $X$ alone.

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.