ccoeft.f

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

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

Variables on this page (5)

a, b, c, d, json

Purpose

This file sets up the tridiagonal system for one implicit time step of composition transport by rotationally induced mixing, on an evenly spaced grid in the coordinate $\chi$ built by GETGRID and MIXGRID (the comments' "radius" and "cm" are older than that grid).

The comments describe an older use. They speak of an even grid in radius, a spacing in cm, coefficients "COMPUTED IN SR COEFFT" (lines 25–36 ↗) and "ANGULAR MOMENTUM" at the boundary (lines 61–63 ↗), but MIXCOM calls the routine with DCHI, the spacing of the grid in $\chi$, for the composition (mixcom.f, lines 84–89 ↗), with coefficients from CODIFF and MIXGRID, so only the comments are out of date.

How it fits

Inputs: the arguments COD (diffusion coefficients, "computed in SR COEFFT" per the comment), DR (grid spacing: MIXCOM passes DCHI, the spacing in $\chi$), DT (time step, s), ECOMP (mass fractions at the start of the step on the even grid), EM (the masses of the grid's shells, g), NTOT (number of points).
Who calls it: MIXCOM.
Outputs: the coefficients A, B, C and right-hand side D in TRIDI_MOD, which CTRIDI then solves (MIXCOM calls only that one); TRIDI_MOD is shared with these and with DCOEFT and GETQUA.

The method

The comment gives the equation as $$\frac{\partial X}{\partial t} = \frac{1}{4\pi\rho r^2}\frac{\partial}{\partial r}\left(D\,\frac{\partial X}{\partial r}\right),$$ differenced fully implicitly on the grid in $\chi$ ($\Delta r$ below is DR, that is, $\Delta\chi$; $m_i$ is the shell mass EM): $$X_i^{n+1} - X_i^{n} = \frac{\Delta t}{\Delta r\,m_i}\left[D_{i+1}\left(X_{i+1}^{n+1} - X_i^{n+1}\right) - D_i\left(X_i^{n+1} - X_{i-1}^{n+1}\right)\right],$$ where COD(I) is $D$ at the interface between points $i-1$ and $i$ (written $D_{i-1/2}$ in the comment). So, with $f_i = \Delta t/(\Delta r\, m_i)$: $$a_i = -f_i D_i, \qquad b_i = 1 + f_i (D_i + D_{i+1}), \qquad c_i = -f_i D_{i+1}, \qquad d_i = X_i^n.$$ At both ends there is no flux across the boundary: $a_1 = 0$, $b_1 = 1 + f_1 D_2$; $c_N = 0$, $b_N = 1 + f_N D_N$.

The scheme conserves mass: CODIFF includes $4\pi$ in COD2, and MIXGRID multiplies the coefficient by $\rho r^2\,d\chi/dr$, so COD is $4\pi\rho r^2 D\,d\chi/dr$, and the flux between shells is COD times the difference in $X$ over $\Delta\chi$.

Definitions

subroutineCCOEFT

lines 2–82 · GitHub

Builds the implicit tridiagonal diffusion system for composition on the even grid in $\chi$, with zero-flux boundaries, into TRIDI_MOD.

Uses modules DIMENSIONS_MOD, TRIDI_MOD

Called by MIXCOM

Change log

The changes since rp_1

Changed since rp_1. Declarations only: /TRIDI/ to its module, the unused /CCOUT/ and /LUOUT/ dropped, IMPLICIT NONE.

1C $$$$$$ 2 SUBROUTINE CCOEFT(COD,DR,DT,ECOMP,EM,NTOT) 3C 10/15/26 MHP via Claude 4C code changed to take JSON from DIMENSIONS_MOD instead of declaring 5C PARAMETER(JSON=5000) here. 6 USE DIMENSIONS_MOD 7C 10/14/26 MHP via Claude 8C code changed to convert COMMON/TRIDI/ into TRIDI_MOD. 9 USE TRIDI_MOD 10C 10/20/26 MHP via Claude 11C code changed to replace implicit typing with IMPLICIT NONE and 12C explicit declarations. 13 IMPLICIT NONE 14 REAL*8 COD,DR,DT,ECOMP,EM,FACT,FACT0 15 INTEGER I,NTOT 16 17C 10/13/26 MHP via Claude: removed the unused COMMON/NMLPAR/ declaration 18C (no member referenced here; residue of the obsolete-option removals). 19C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration 20C (no member referenced here; residue of the obsolete-option removals). 21C 10/13/26 MHP via Claude: removed the unused COMMON/IOUNITS/ declaration 22C (no member referenced here; residue of the obsolete-option removals). 23 DIMENSION COD(JSON),ECOMP(JSON),EM(JSON) 24 SAVE 25C SET UP COEFFICIENTS OF DIFFUSION DIFFERENCE EQUATIONS FOR COMPOSITION 26C TRANSPORT DUE TO ROTATIONALLY INDUCED MIXING. 27C 28C INPUT VARIABLES: 29C COD : RUN OF DIFFUSION COEFFICIENTS (COMPUTED IN SR COEFFT). 30C DR : GRID SPACING IN RADIUS(CM). 31C DT : DIFFUSION TIMESTEP (SEC). 32C ECOMP : STARTING ARRAY OF COMPOSITION MASS FRACTION AT THE EQUALLY 33C SPACED GRID. 34C EM : RUN OF MASSES OF THE EQUALLY SPACED GRID POINTS (GM). 35C NTOT : NUMBER OF POINTS IN THE EQUALLY SPACED GRID USED TO SOLVE 36C THE DIFFUSION EQUATION. 37C OUTPUT VARIABLES : 38C THE ABUNDANCE OF SHELL I AT TIME N+1 (COMP(I,N+1)) IS A FUNCTION OF 39C COMP(I-1),COMP(I),AND COMP(I+1) AS DISCUSSED BELOW. THIS CAN BE EXPRESSED 40C AS A TRIDIAGONAL MATRIX. TERMS WHICH DEPEND ON COMP(I-1) ARE STORED IN 41C ARRAY A; THOSE WHICH DEPEND ON COMP(I) ARE STORED IN ARRAY B; THOSE 42C WHICH DEPEND ON COMP(I+1) ARE STORED IN ARRAY C; AND THE STARTING RUN 43C OF COMPOSITION IS STORED IN ARRAY D. 44C 45C THE DIFFUSION EQUATION WE ARE SOLVING IS 46C dX/dT = 1/(4pi*RHO*R**2)d/dR(D*dX/dR) 47C WHERE D IS OUR DIFFUSION COEFFICIENT,X IS THE MASS FRACTION OF THE 48C SPECIES IN QUESTION,T IS TIME,R IS RADIUS, AND RHO IS DENSITY. 49C USING M = 4pi*RHO*R**2*DR THIS IS DIFFERENCED FOR SHELL I AS 50C X(I,N+1)-X(I,N)=(DT/DR)*(1/M)*(D(I+1/2)*(X(I+1,N+1)-X(I,N+1)) - 51C D(I-1/2)*(X(I,N+1)-X(I-1,N+1))) 52C WHERE M IS THE MASS OF SHELL , INDEX N+1 REFERS TO VALUES AT THE 53C END OF THE TIMESTEP, INDEX N REFERS TO VALUES AT THE BEGINNING OF 54C THE TIMESTEP. 55C THIS SYSTEM IS SUPPLEMENTED WITH APPROPRIATE BOUNDARY CONDITIONS. 56C 57C STORE THE RUN OF COMPOSITION AT THE BEGINNING OF THE TIMESTEP IN ARRAY D 58 DO 10 I = 1,NTOT 59 D(I) = ECOMP(I) 60 10 10 CONTINUE 61C FIRST SHELL BOUNDARY CONDITIONS: NO FLOW BELOW THE BOUNDARY 62C I.E. THE ANGULAR MOMENTUM TRANSPORT AT THE FIRST SHELL DOES NOT 63C DEPEND ON THE SHELLS BELOW IT.(SAME USED FOR LAST SHELL B.C.) 64 FACT0 = DT/DR 65 FACT = FACT0/EM(1) 66 A(1) = 0.0D0 67 B(1) = 1.0D0 + FACT*COD(2) 68 C(1) = -FACT*COD(2) 69 DO 20 I = 2,NTOT-1 70 FACT = FACT0/EM(I) 71 A(I) = -FACT*COD(I) 72 B(I) = 1.0D0 + FACT*(COD(I)+COD(I+1)) 73 C(I) = -FACT*COD(I+1) 74 20 20 CONTINUE 75C LAST SHELL BOUNDARY CONDITIONS. 76 FACT = FACT0/EM(NTOT) 77 A(NTOT) = -FACT*COD(NTOT) 78 B(NTOT) = 1.0D0 + FACT*COD(NTOT) 79 C(NTOT) = 0.0D0 80 81 RETURN 82 END

1"""Python version of ccoeft.f (YREC 7be4133): the implicit tridiagonal 2system for one step of composition mixing on the even grid.""" 3 4from pyrec import tridi_mod 5import numpy as np 6 7 8def ccoeft(cod: np.ndarray, dr: float, dt: float, ecomp: np.ndarray, em: np.ndarray, ntot: int) -> None: 9 """Fill the tridiagonal system a, b, c, d of tridi_mod for one fully 10 implicit diffusion step of a mass fraction on the even grid in chi, 11 with no flux across either end. 12 13 Translates: ccoeft.f#CCOEFT (7be4133) 14 15 Arguments: cod, the diffusion coefficients at the interfaces (cod[i] 16 between points i-1 and i); dr, the grid spacing (MIXCOM passes the 17 spacing in chi); dt, the time step (s); ecomp, the mass fractions at 18 the start of the step; em, the shell masses (g); ntot, the number of 19 points. None is set; the results go to tridi_mod. 20 21 """ 22 23 a = tridi_mod.a 24 b = tridi_mod.b 25 c = tridi_mod.c 26 d = tridi_mod.d 27 28 # The composition at the start of the step (Fortran lines 58–60). 29 for i in range(1, ntot + 1): 30 d[i] = ecomp[i] 31 32 # First shell: no flow below it (Fortran lines 64–68). 33 fact0 = dt / dr 34 fact = fact0 / em[1] 35 a[1] = 0.0 36 b[1] = 1.0 + fact * cod[2] 37 c[1] = -fact * cod[2] 38 39 # The interior shells (Fortran lines 69–74). 40 for i in range(2, ntot): 41 fact = fact0 / em[i] 42 a[i] = -fact * cod[i] 43 b[i] = 1.0 + fact * (cod[i] + cod[i + 1]) 44 c[i] = -fact * cod[i + 1] 45 46 # Last shell: no flow above it (Fortran lines 76–79). 47 fact = fact0 / em[ntot] 48 a[ntot] = -fact * cod[ntot] 49 b[ntot] = 1.0 + fact * cod[ntot] 50 c[ntot] = 0.0