fix low density/Pressure SCVH EOS tables

[Debraheem]

This pr is an attempt to address #995, which spawned from a mailing-list thread initially raised by Natasha Ivanova:

I believe the main concern was the mean molecular weight being that of an ionized mixture in the regions highighted in the plots below, when adopting both Free_EOS and OPAL/SCVH.

https://lists.mesastar.org/pipermail/mesa-users/2026-May/015873.html:

Particularly in the low T-rho regime, as I'm not sure there is a better fallback than HELM in the high pressure regime (that would require an actual eos for that region).

For the SCVH tables, below logP_min = -0.6 the SCVH coverage is now supplied analytically in
eos/eosDT_builder/src/scvh_core.f90. For each fixed logT row and each species, the tail uses the edge state and extrapolates to the ideal gas limit.

logP0 = -0.6
logrho0 = logrho(logP0, logT)
logu0 = logu(logP0, logT)
logs0 = logs(logP0, logT)
R = 10**(logP0 - logrho0 - logT)

and the ideal low-pressure limit

logrho_ideal(logP) = logrho0 + logP - logP0
logu_ideal(logP) = logu0
S_ideal(logP) = S0 + R*ln(10)*(logP0 - logP)
logs_ideal(logP) = log10(S_ideal)
comp_ideal(logP) = comp0

The implementation uses this ideal value below the first lower-pressure grid
step and applies a local cubic correction between logP0 and that point so the
tail is value-continuous and has the same d/dlogP as the restored processed
SCVH table at the edge. That keeps the generated eosDT interpolation support
smooth without editing the raw SCVH data files.

In this pr, I also noticed when regenerating the eos tables that there was a HELM fallback in
the high pressure region, and it produced noisy mu features
near the OPAL/SCVH support boundary in regenerated tables. I think this fallback didn't exist,
when the tables where generated the first time, so i turned it off for the table generation. see check_results :

was
gamma1_opalscvh <= 1d0 .or. gamma1_opalscvh >= 2d0 .or. grad_ad_opalscvh >= 50d0) then
now
gamma1_opalscvh <= 0d0) then

This might not be correct, but it was the only way for me to make the new eos smooth like the old eos tables, but maybe this expression should be adjusted.

See the plots below

After adding the smooth extrapolation to SCVH, i also noticed there was a bad row of partials on the
FreeEOS SCVH boundary at logT = 3 - 3.1. It turns out the fallback here were rows generated from the
MESA EOS with lnPgas, lnE, and lnS in natural log form, but the FreeEOS data files
store those columns as base-10 logarithms. so i these just need to be converted in eos/eosFreeEOS_builder/src/free_eos_table.f90. (See the original issue below).

I had to regenerate the eos for this pr.

• eos/eosFreeEOS_data.tar.xz
• eos/eosDT_data.tar.xz

I bumped the the FreeEOS table version and the eosDT table version to 52.

A few plots from the new eos looks like this:

(a bit low res for the eos plotter)

The eosDT region definitions shifted very slightly in eos/eosDT_builder/src/eos_regions_defs.dek and
eos/eosDT_builder/src/eos_regions_code.dek keep the generated OPAL/SCVH support
out of the warm, very-low-density FreeEOS side. This prevents the SCVH
low-pressure floor from interacting with the FreeEOS boundary at
logrho ~ -15.

This PR probably needs some further attention validation and testing, maybe clean up.

[evbauer]

Wow. Did you actually rerun the full FreeEOS preprocessor to regenerate all of the FreeEOS tables after the fix for the MESA fallback points? What version of FreeEOS did you use? Last time I tried this, I was having trouble getting an external FreeEOS installation to compile, run, and be compatible with the MESA preprocessor. But that was a while ago.

[evbauer]

Rather than extrapolating "OPAL/SCVH" over to low density, why not just blend over to neutral ideal gas there? It seems like making a table that extrapolates toward ideal gas is maybe just a less-optimal version, since we can do it analytically with auto-diff now?

[evbauer]

Great catch! I think you're right about this. Do you know if it has an impact on anything?

In practice, I think we cut out most of the region where FreeEOS was failing and falling back to MESA with the "logQ_cut"-related boundary controls. Anyway, I've always been a bit uncomfortable with the idea that our "FreeEOS" tables that were meant to replace the older EOS might end up falling back to the older EOS without being very clear about that. So I like the fix, but I don't really like the existence of the thing that we're fixing here...

[Debraheem]

I've definitely seen models crash from this, in fact this issue from jaime comes to mind : #822

[Debraheem]

after asking an ai to check this further:

It affects every row marked MESA = 1 in the FreeEOS tables, since those rows were written through the mesa_eos_eval fallback path. It only showed up at the low-T boundary because that was the place where runtime interpolation of the FreeEOS data was actually reaching those MESA = 1 rows.

total numeric rows: 14,399,221
MESA fallback rows: 2,618,641

examples pulled from the tables comparing the old and new free eos tables:

Here are a few matched v51/v52 fallback rows from mesa-FreeEOS_02z70x.data. In all of these rows, the v51 values are larger than the v52 values by ln(10), which is the sign that v51 wrote natural-log values into fields that are read as base-10 logs.

logT	T [K]	logRho	rho [g/cm^3]	v51 logP/logE/logS	v52 logP/logE/logS
2.10	1.26e2	-17.89	1.29e-18	-10.0294 / 22.6406 / 20.9587	-4.3557 / 9.8327 / 9.1022
2.50	3.16e2	-9.20	6.31e-10	1.9463 / 23.8743 / 20.7341	0.8453 / 10.3685 / 9.0047
3.00	1.00e3	-0.37	4.27e-1	23.4272 / 24.6846 / 19.4259	10.1743 / 10.7204 / 8.4366
3.50	3.16e3	0.18	1.51e0	29.2003 / 28.4440 / 19.7273	12.6815 / 12.3531 / 8.5675
4.00	1.00e4	1.09	1.23e1	33.3487 / 30.8600 / 19.6948	14.4831 / 13.4023 / 8.5533
5.00	1.00e5	1.65	4.47e1	35.7939 / 32.1696 / 20.1509	15.5451 / 13.9711 / 8.7514

|

[evbauer]

Why did the write format change?

[Debraheem]

Some values overflowed the old fixed-width f10.5 columns when I regenerated the eosDT tables as part of this PR. I still need to check why this overflow happened. It might be symptoms of another bug i might have introduced in this pr, so i'm looking closer.

[Debraheem]

Wow. Did you actually rerun the full FreeEOS preprocessor to regenerate all of the FreeEOS tables after the fix for the MESA fallback points? What version of FreeEOS did you use? Last time I tried this, I was having trouble getting an external FreeEOS installation to compile, run, and be compatible with the MESA preprocessor. But that was a while ago.

No, I did not rerun the external FreeEOS code or regenerate the raw FreeEOS data. The raw FreeEOS part is the same. The confusing bit is that MESA’s eosFreeEOS data is not purely raw FreeEOS output. It also contains rows where the FreeEOS builder and falls back to the MESA EOS, marked by MESA = 1. This PR patches/re-stitches those fallback rows using the corrected MESA EOS behavior, including the SCVH low pressure extension, and fixes the log base for lnPgas, lnE, and lnS in those fallback rows. So this is not a new raw FreeEOS table generation. It is a correction to the MESA-packaged FreeEOS support tables in the regions where they were already using MESA fallback values.

I have a longer term goal of replacing this entire region with pure FreeEOS using FreeEOS3 and maybe chemEOS at lowT in a Free Energy formulation (your initial idea @evbauer), but I can't do this at the moment. Indeed, I also had the same issues as you regenerating the original tables last i tried, but I think I've learned enough now to give it another attempt soon. Perhaps with AI and some of @VincentVanlaer's compiler knowledge we can get FreeEOS to build raw tables again in the future.

Rather than extrapolating "OPAL/SCVH" over to low density, why not just blend over to neutral ideal gas there? It seems like making a table that extrapolates toward ideal gas is maybe just a less-optimal version, since we can do it analytically with auto-diff now?

This is close to what I did, but not quite the same as blending to the existing neutral ideal EOS at runtime.

In this PR the low-P SCVH extension is built during eosDT preprocessing. Below the processed SCVH edge at logP_min = -0.6, the code makes an edge-matched ideal-gas tail and writes that into the generated OPAL/SCVH support tables. The thermodynamic variables are pushed toward the ideal low-pressure behavior, but the H/He molecular and ionization fractions are taken from the SCVH table edge and then held fixed in pressure through the tail. See scvh_core.f90 around the new interp_vals_lowP_tail/tail_species logic, especially the calls that carry compv1/compv2 from the edge, and the later mu/logNe calculation.

So this is not exactly a pure neutral ideal blend. A runtime blend to the existing ideal EOS would be cleaner architecturally, but it would also change the EOS blending machinery and would give neutral-ideal composition unless we added a separate analytic treatment for H2/H/H+ and He/He+/He++. I did it in the preprocessor because it is a more local fix: the generated eosDT table stays smooth at the old SCVH edge, the region stops falling back to HELM, and the runtime EOS selection does not need a new component boundary.

on why i did it in the tables as opposed to doing it at runtime during the eos evaluations:

My fix did not need to be in the tables. The same tail formula could be evaluated analytically. However, in MESA’s current architecture, OPAL/SCVH at runtime is already reading in a preprocessed eosDT table. So putting the tail in the preprocessor is the lowest-impact place to fix the existing generated table and not introducing a bunch of other heuristics that have to computed at runtime.

The distinction is:

• Your suggestion, an analytic runtime version: add a new runtime branch/blend that detects the low-density SCVH edge and evaluates an analytic tail or neutral ideal EOS there.

• My current PR: do the analogous extension in the preprocessing stage, so the generated OPAL/SCVH eosDT tables have an analytic edge-matched tail. Runtime sees a smooth table and does not need new EOS selection logic.

The runtime version is possible, but it's a bigger change. If it used the existing neutral ideal EOS, it would not reproduce what I am doing either, since the current PR is keeping the ionization fractions anchored to the SCVH edge state (H2/H/H+ and He/He+/He++ fractions). The analytic approach at evaluation time could grab the SCVH edge state as a function of logT, then use the requested logP to carry those fractions into the tail and match energy/entropy/derivatives at the edge every time the routine is called. But then you are effectively moving the preprocessing logic into runtime, and recomputing that matching every call. At that point it seemed simpler to put the smooth extension into the table once and let the existing eosDT machinery handle it. And then there is no confusion in scenarios when someone looks at the Free_EOS data file in the data directory they arent pulling out Helm values and getting confused, because the precprocessor was falling back to HELM (this is what's happening now).

So the pragmatic reason for tables is: the existing OPAL/SCVH product is a table, and the problem was just confusing table support/fallback behavior that's not communicated to the user (so the free eos data files have a lot of fallback from MESA's other EOS'). Replacing the HELM fallback smoothly in preprocessing fixes the immediate issue without inventing a new runtime EOS boundary.

Maybe your analytic runtime approach could be cleaner, but i figured this pr is smaller in what it actually changes. I'm happy with whatever is ultimately cleaner/better. I don't mind backing up and switching the approach. Also maybe it's not necessary to keep the ionization states, in which case an ideal eos blend might be better? What do you think, i defer to your judgement on the best path forward.

[Debraheem]

I had some AI help going through the differences, and I made plots comparing the tables in this PR to the tables currently on main. One useful detail: in git history, eos/eosDT_data.tar.xz is unchanged between the initial MESA r15140 import and this PR, so I cannot tell from git when the v51 tables were originally generated. They predate the public GitHub history.

The main difference I found is in the eosDT preprocessor path. The v51/main generation source uses:

scvh_only = .true.

so the builder forces SCVH first across the eosDT output grid. It then applies the stricter check_results filter, including:

gamma1_opalscvh <= 1d0 .or. gamma1_opalscvh >= 2d0 .or. grad_ad_opalscvh >= 50d0

If that forced-SCVH point fails the filter, the driver retries with OPAL; if that also fails, it falls back to HELM.

When I regenerated the tables, the OPAL/SCVH region looked noisier than the v51 tables in this region:

The fix that made to make the regenerated tables look smooth again was to stop forcing SCVH everywhere and instead use:

opal_scvh_only = .true.
scvh_only = .false.

so opal_scvh_driver can use its internal OPAL/SCVH boundary/blend. I also relaxed the gamma/grad_ad part of check_results from the strict condition above to:

gamma1_opalscvh <= 0d0

while still rejecting bad values. I feel like relaxing this only made the gamma1 plots smoother on the new branch as compared to main.

So I think the difference right is:

v51/main behavior: force SCVH, reject suspicious gamma1/grad_ad values, patch rejected points with OPAL
PR behavior: use normal OPAL/SCVH selection/blending, reject only clearly broken OPAL/SCVH points

The high-density differences and f10.5 overflow are a separate side effect of changing that preprocessor path. Some fallback points that were probably OPAL-patched in the original v51-preprocessor path can become HELM fallback points in the new path. Current HELM returns finite but very large derivatives at a few of those high-density support points, which exposed the fixed width write format problem. THat occurs in the

The derivative plots are sensitive to all of this because MESA builds bicubic interpolation coefficients from the regenerated table surface. So changing whether a patch is SCVH, OPAL, OPAL/SCVH blend, HELM, or the new low-P tail can change d/dlnT and d/dlnRho even when the value plots look broadly similar.

Here are some comparison plots to main. Somethings are improved, others maybe a little worse. curious for feedback.
compare_branch_to_main.pdf

Either way, the main lesson im learning is that the eos data tables haven't been regenerated in a long time, and some things have changed since then. I'm slowly going through all the noisy patches to see if i can quench those as well.