rhoofp06.f

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

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

Variables on this page (10)

eos, mv, mx, nr, nra, nt, rho, t6list, xa, xz

Purpose

Inverts the OPAL 2006 equation of state: the density for a given $X$, $T_6$ and pressure. It follows the method of RHOOFP and RHOOFP01, the versions for the original and 2001 tables (both removed since rp_1), with a different treatment of radiation.

How it fits

Called by OEQOS06; calls ESAC06 for each trial density. Modules: AEOS06_MOD (tables and grids from READCOEOS06) and EOS06MISC_MOD (EOS, the results of ESAC06, and NRA, the last density index at each temperature, set by a DATA statement in the module and also used by EQBOUND06). It returns the density, or $-999$.

The method

Find the bracketing composition index MLO in XA and temperature index KLO in T6LIST(1,·), by bisection.
Estimate the table's pressure range at that temperature from the first stored quantity, read as $P/(\rho T_6)$: $P_\mathrm{max} = \texttt{XZ}(\cdot,1,k,n_k)\,T_6\,\rho_{n_k}$ and $P_\mathrm{min} = \texttt{XZ}(\cdot,1,k,1)\,T_6\,\rho_1$, with $n_k$ = NRA(KLO). If $P_\mathrm{gas} > 1.25\,P_\mathrm{max}$ or $P_\mathrm{gas} < P_\mathrm{min}$, it fails.
First guess $\rho_1 = \rho_{n_k} P_\mathrm{gas}/P_\mathrm{max}$; a second guess five times larger (capped at $\rho_{n_k}$) or five times smaller (floored at $10^{-14}$), depending on the side.
Up to 11 steps of the secant method, keeping a bracket (regula falsi): $$\rho_3 = \rho_1 + (\rho_2 - \rho_1)\,\frac{P_\mathrm{gas} - P_1}{P_2 - P_1}.$$

Failures (out of the table, no convergence, an error return from ESAC06) all return $-999$ silently; the messages are commented out.

Differences from the removed RHOOFP.

Radiation. When irad is 1, the radiation pressure is subtracted first, $$P_\mathrm{gas} = P - \tfrac43\,\frac{\sigma}{c}\,T_6^4\quad(\text{in Mb}),$$ and the search matches $P_\mathrm{gas}$ against ESAC06 called with irad = 0 (gas only), with $P_\mathrm{min}$ and $P_\mathrm{max}$ also without radiation. The other two versions matched the total pressure instead.
The edge NRA comes from EOS06MISC_MOD, not from a local DATA list.
Tolerance. Convergence is $|P_3 - P|/P < 10^{-5}$ (lines 117–121 ↗), changed from $5\times10^{-8}$ "to prevent array value eos(5) from growing without bound and crashing the program during certain model runs" (MR 2025-10-10). The growth comes from ESAC06: called here with IORDER = 1 and IRAD = 0, it rescales EOS(5), which it has not interpolated, on every call, so fewer calls mean less growth (see esac06.f, suspected bug 1).

Question 1 (line 120 ↗). The looser tolerance. The 2025 change accepts a density whose pressure is off by up to $10^{-5}$ relative. Is that accurate enough for the hydrostatic solution, and should the stricter tolerance come back once ESAC06 scales only what it interpolates?

Definitions

functionRHOOFP06

lines 8–147 · GitHub

Uses modules AEOS06_MOD, EOS06MISC_MOD

Called by OEQOS06

Flowchart of its calls, in order (4)

ESAC06 91

if (p1 .gt. pnr), lines 93–100

ESAC06 99

else, lines 101–108

ESAC06 106

ESAC06 115

Compiler warnings (3)

gfortran -Wall -Wextra, compiling only; nothing was run.

line 20: Unused parameter ‘mv’ declared [-Wunused-parameter]
line 20: Unused parameter ‘nr’ declared [-Wunused-parameter]
line 66: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. The first-call switch (ITIME in /LREADCO/) and the call to ESAC06 at a fixed point that forced the tables to be read are removed (the tables are now read once by SETEOS06, from SETUPS); SIGMACC is a double-precision constant (1.8914785D-3), which settles former question 3; RHOOFP and RHOOFP01 are removed, which settles former question 2 (making the three versions agree); COMMON blocks to modules, IMPLICIT NONE.

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C 3C RHOOFP06 4C 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 7C function rhoofp06(x,ztab,t6,p,irad) ! KC 2025-05-31 8 function rhoofp06(x,t6,p,irad) 9C 10/20/26 MHP via Claude 10C code changed to convert COMMON/AEOS06/, COMMON/EOS06MISC/ into 11C AEOS06_MOD, EOS06MISC_MOD. 12 USE AEOS06_MOD 13 USE EOS06MISC_MOD 14C 10/20/26 MHP via Claude 15C code changed to replace implicit typing with IMPLICIT NONE and 16C explicit declarations. 17 IMPLICIT NONE 18 REAL*8 P,P1,P2,P3,PMAX,PMIN,PNR,PR,RAT,RHOG1,RHOG2,RHOG3,RHOOFP06, 19 * * SIGMACC,T6,X 20 INTEGER ICOUNT,IHI,ILO,IMD,IRAD,KLO,MLO,MV,MX,NR,NT 21 22 parameter (mx=5,mv=10,nr=169,nt=197) 23C 10/13/26 MHP via Claude 24C code changed to remove COMMON/LREADCO/; the table is now read once by 25C SETEOS06 from SETUPS, so this routine no longer needs the first-time 26C switch or the fabricated ESAC06 call that used to force the read. 27C 10/9/26 MHP via Claude 28C code changed to rename AEOS06's M member to MIDX: M is such a 29C universally-used name elsewhere in this codebase (e.g. shell count 30C in calcad.f) that keeping it here was just asking for a future 31C collision. M is genuinely used (as a composition/mixture-axis 32C index) in esac06.f and readcoeos06.f; declared but never read or 33C written by name in this file. 34C KC 2025-05-30 reordered common block elements 35C common/beos06/ iri(10),index(10),nta(nr),zz(mx),nra(nt) 36C 10/9/26 MHP via Claude 37C code changed to merge COMMON/EEOS06/, COMMON/AAEOS06/, COMMON/BBEOS06/, 38C and COMMON/BEOS06/ into COMMON/EOS06MISC/ (see esac06.f for the full 39C rationale). AEOS06 (above) is left untouched. 40 data sigmacc/1.8914785D-3/ 41c-------------------------------------------------------------------- 42 43 save 44 rat=sigmacc 45 pr=0D0 46 if(irad .eq. 1) pr=4D0/3D0*rat*t6**4 ! Mb 47 pnr=p-pr 48 49 ilo=2 50 ihi=mx 51 8 8 if(ihi-ilo .gt. 1) then 52 imd=(ihi+ilo)/2 53 if(x .le. xa(imd)+1.D-7) then 54 ihi=imd 55 else 56 ilo=imd 57 endif 58 go to 8 59 endif 60 mlo=ilo 61 62 ilo=nt 63 ihi=2 64 11 11 if(ilo-ihi .gt. 1) then 65 imd=(ihi+ilo)/2 66 if (t6 .eq. t6list(1,imd)) then 67 ilo=imd 68 go to 14 69 endif 70 if(t6 .le. t6list(1,imd)) then 71 ihi=imd 72 else 73 ilo=imd 74 endif 75 go to 11 76 endif 77 14 14 klo=ilo 78 79 pmax=xz(mlo,1,klo,nra(klo))*t6*rho(nra(klo)) 80 pmin=xz(mlo,1,klo,1)*t6*rho(1) 81 if((pnr .gt. 1.25D0*pmax) .or. (pnr .lt. pmin)) then 82c write (ISHORT,'(" The requested pressure-temperature not in", 83c x " the OPAL 2006 EOS table")') 84c stop 85c write (ISHORT,'("pnr, pmax,pmin=",3e14.4)') pnr,pmax,pmin 86 goto 999 !RHOOFP06 error exit 87 endif 88 89 rhog1=rho(nra(klo))*pnr/pmax 90C call esac06 (x,ztab,t6,rhog1,1,0,*999) ! KC 2025-05-31 91 call esac06 (x,t6,rhog1,1,0,*999) 92 p1=eos(1) 93 if(p1 .gt. pnr) then 94 p2=p1 95 rhog2=rhog1 96 rhog1=0.2D0*rhog1 97 if(rhog1 .lt. 1.D-14) rhog1=1.D-14 98C call esac06 (x,ztab,t6,rhog1,1,0,*999) ! KC 2025-05-31 99 call esac06 (x,t6,rhog1,1,0,*999) 100 p1=eos(1) 101 else 102 rhog2=5D0*rhog1 103C if(rhog2 .gt. rho(klo)) rhog2=rho(klo) ! Corrected below llp 8/19/08 104 if(rhog2 .gt. rho(nra(klo))) rhog2=rho(nra(klo)) ! Had wrong pointer, see rhog1= ten lines up 105C call esac06 (x,ztab,t6,rhog2,1,0,*999) ! KC 2025-05-31 106 call esac06 (x,t6,rhog2,1,0,*999) 107 p2=eos(1) 108 endif 109 110 icount=0 111 1 1 continue 112 icount=icount+1 113 rhog3=rhog1+(rhog2-rhog1)*(pnr-p1)/(p2-p1) ! KC 2025-05-31 114C call esac06 (x,ztab,t6,rhog3,1,0,*999) 115 call esac06 (x,t6,rhog3,1,0,*999) 116 p3=eos(1) 117C Changed the comparison below to use the commented-out value 1.D-5 118C found here to prevent array value eos(5) from growing without bound 119C and crashing the program during certain model runs. - MR 2025-10-10 120 if (abs((p3-pnr)/pnr) .lt. 1.D-5) then 121C IF (DABS((P3-PNR)/PNR) .LT. 0.5D-7) THEN 122 rhoofp06=rhog3 123 return 124 endif 125 if (p3 .gt. pnr) then 126 rhog2=rhog3 127 p2=p3 128 if (icount .lt. 11) go to 1 129c write (ISHORT,'("Rhoofp06: No convergence after 10 tries")') 130 goto 999 !RHOOFP06 error exit 131c stop 132 else 133 rhog1=rhog3 134 p1=p3 135 if (icount .lt. 11) go to 1 136c write (ISHORT,'("RHOOFP06: No convergence after 10 tries")') 137 goto 999 ! RHOOFP06 error exit 138c stop 139 endif 140 141 999 999 CONTINUE 142 RHOOFP06=-999.0D0 143C WRITE(ISHORT,'("RHOOFP06: FAILED TO FIND RHO")') 144 RETURN 145 146 147 end 148

1"""Python version of rhoofp06.f (YREC 7be4133): the density for a given X, 2T6 and pressure, from the OPAL 2006 equation of state.""" 3 4from pyrec import aeos06_mod 5from pyrec import eos06misc_mod 6from pyrec.esac06 import esac06 7 8MX = 5 9MV = 10 10NR = 169 11NT = 197 12SIGMACC = 1.8914785e-3 # DATA: a/3 in the table's units (Mb per T6**4) 13 14 15def rhoofp06(x: float, t6: float, p: float, irad: int) -> float: 16 """Invert the OPAL 2006 EOS: find the density whose gas pressure 17 matches p (less the radiation pressure when irad is 1), by bisection 18 in the table's X and T6 grids, then up to 11 steps of regula falsi on 19 esac06 called without radiation. 20 21 Translates: rhoofp06.f#RHOOFP06 (7be4133) 22 23 Arguments: x (X), t6 (T in 1e6 K), p (pressure, Mb), irad (1: p 24 includes radiation). None is set. 25 26 Returns the density, or -999.0 on failure (out of the table, no 27 convergence, or esac06's alternate return). 28 29 """ 30 31 # The gas pressure (Fortran lines 44–47). 32 a = aeos06_mod 33 eos = eos06misc_mod.eos 34 nra = eos06misc_mod.nra 35 rat = SIGMACC 36 pr = 0.0 37 if irad == 1: 38 pr = 4.0 / 3.0 * rat * t6**4 # Mb 39 pnr = p - pr 40 41 # Bisection for the composition index (label 8) (Fortran lines 49–60). 42 ilo = 2 43 ihi = MX 44 while ihi - ilo > 1: 45 imd = (ihi + ilo) // 2 46 if x <= a.xa[imd] + 1.0e-7: 47 ihi = imd 48 else: 49 ilo = imd 50 mlo = ilo 51 52 # Bisection for the temperature index (label 11), from the top (Fortran 53 # lines 62–77). 54 ilo = NT 55 ihi = 2 56 while ilo - ihi > 1: 57 imd = (ihi + ilo) // 2 58 if t6 == a.t6list[1, imd]: 59 ilo = imd 60 break # go to 14 61 if t6 <= a.t6list[1, imd]: 62 ihi = imd 63 else: 64 ilo = imd 65 klo = ilo # label 14 66 67 # The table's pressure range at this temperature (Fortran lines 79–87). 68 pmax = a.xz[mlo, 1, klo, nra[klo]] * t6 * a.rho[nra[klo]] 69 pmin = a.xz[mlo, 1, klo, 1] * t6 * a.rho[1] 70 if pnr > 1.25 * pmax or pnr < pmin: 71 return -999.0 # go to 999 72 73 # Two first guesses (Fortran lines 89–108). 74 rhog1 = a.rho[nra[klo]] * pnr / pmax 75 if esac06(x, t6, rhog1, 1, 0): 76 return -999.0 77 p1 = eos[1] 78 if p1 > pnr: 79 p2 = p1 80 rhog2 = rhog1 81 rhog1 = 0.2 * rhog1 82 if rhog1 < 1.0e-14: 83 rhog1 = 1.0e-14 84 if esac06(x, t6, rhog1, 1, 0): 85 return -999.0 86 p1 = eos[1] 87 else: 88 rhog2 = 5.0 * rhog1 89 if rhog2 > a.rho[nra[klo]]: 90 rhog2 = a.rho[nra[klo]] 91 if esac06(x, t6, rhog2, 1, 0): 92 return -999.0 93 p2 = eos[1] 94 95 # Regula falsi (label 1), at most 11 steps (Fortran lines 110–139). 96 icount = 0 97 while True: 98 icount = icount + 1 99 rhog3 = rhog1 + (rhog2 - rhog1) * (pnr - p1) / (p2 - p1) 100 if esac06(x, t6, rhog3, 1, 0): 101 return -999.0 102 p3 = eos[1] 103 # Question (see the page): the tolerance 1e-5 (MR 2025-10-10). 104 if abs((p3 - pnr) / pnr) < 1.0e-5: 105 return rhog3 106 if p3 > pnr: 107 rhog2 = rhog3 108 p2 = p3 109 else: 110 rhog1 = rhog3 111 p1 = p3 112 if icount >= 11: 113 return -999.0 # no convergence