fix low density/Pressure SCVH EOS tables

atm/test/test_output

@@ -85,17 +85,17 @@
 
  test_irradiated
 
-                              test_grey_irradiated: kap    4.3649748642304882D-03
+                              test_grey_irradiated: kap    4.3649748858080083D-03
                                             M/M_jupiter    1.4999999999999998D+00
                                             R/R_jupiter    1.5743735095633800D+00
                                                  R/Rsun    1.6178684913857289D-01
                                                   kap_v    4.0000000000000001D-03
                                                       P    1.0000000000000000D+06
 
-                                                    tau    2.9101126947637423D+00
+                                                    tau    2.9101127094055270D+00
                                                    Teff    5.7759999999999891D+03
                                                    T_eq    1.0000000000000000D+03
                                                   T_int    9.0000000000000000D+02
-                                                      T    1.2929349604416070D+03
-                                                   logT    3.1115766787614829D+00
+                                                      T    1.2929349615774468D+03
+                                                   logT    3.1115766791430097D+00
 

docs/source/changelog.rst

@@ -39,6 +39,16 @@ The parameter ``report_max_infall_inside_fe_core`` was ignored in versions r25.1
 
 ``fe_core_infall_limit`` now obeys ``when_to_stop_rtol`` and ``when_to_stop_atol`` again (broken since r11532).
 
+Fixed a bug in the FreeEOS table builder where MESA fallback rows could write
+``lnPgas``, ``lnE``, and ``lnS`` values into table columns that are read as
+base-10 logarithms. This could create discontinuities and noisy EOS partial
+derivatives.
+
+The regenerated eosDT tables now also include an analytic low-pressure SCVH
+tail matched to the SCVH table edge where OPAL/SCVH previously hit table edge
+limits and fell back to HELM. The SCVH region has been extended toward
+:math:`\log_{10}\rho \simeq -14.9`.
+
 
 .. note:: Before releasing a new version of MESA, move `Changes in main` to a new section below with the version number as the title, and add a new `Changes in main` section at the top of the file (see ``changelog_template.rst``).
 

docs/source/known_bugs.rst

@@ -13,6 +13,17 @@ issue, but it may not be complete.
 r26.4.1
 =======
 
+.. _freeeos_fallback_log_units_bug:
+
+EOS: FreeEOS fallback rows use the wrong logarithm base
+------------------------------------------------------
+
+Fixed an issue in the FreeEOS table builder where fallback rows marked by
+``MESA = 1`` could write the MESA EOS values ``lnPgas``, ``lnE``, and ``lnS``
+into table columns that are read as base-10 logarithms. This could create
+discontinuities in the FreeEOS support table and noisy EOS partial derivatives,
+including near the low-density ``logT = 3.0 - 3.1`` FreeEOS boundary.
+
 .. _report_max_infall_inside_fe_core_bug:
 
 Controls: ``report_max_infall_inside_fe_core`` is ignored

eos/eosDT_builder/Makefile

@@ -3,13 +3,13 @@ include ../../make/defaults-module.mk
 # Build
 
 MODULE_NAME := eos-dt-builder
-SRCS := src/azbar.f \
-        src/create_eos_files.f \
+SRCS := src/azbar.f90 \
+        src/create_eos_files.f90 \
         src/eos5_xtrin_h-he.f \
-        src/helm_opal_scvh_driver.f \
+        src/helm_opal_scvh_driver.f90 \
         src/opal_core.f \
-        src/opal_scvh_driver.f \
-        src/scvh_core.f
+        src/opal_scvh_driver.f90 \
+        src/scvh_core.f90
 SRCS_GENERATED := \
         private/helm_alloc.f90 \
         private/helm_polynomials.f90 \
@@ -32,14 +32,14 @@ include ../../make/Makefile
 # Generated sources
 # This just copies over private implementations from eos
 
-$(BUILD_DIR_MODULE)/private/helm_alloc.f90: ../private/helm_alloc.f90 | $(BUILD_DIR_MODULE)/private/
+$(BUILD_DIR_MODULE)/private/helm_alloc.f90: ../private/helm_alloc.f90 | $(BUILD_DIR_MODULE)/private/.
 	@cp $^ $@
 
-$(BUILD_DIR_MODULE)/private/helm_polynomials.f90: ../private/helm_polynomials.f90 | $(BUILD_DIR_MODULE)/private/
+$(BUILD_DIR_MODULE)/private/helm_polynomials.f90: ../private/helm_polynomials.f90 | $(BUILD_DIR_MODULE)/private/.
 	@cp $^ $@
 
-$(BUILD_DIR_MODULE)/private/gauss_fermi.f90: ../private/gauss_fermi.f90 | $(BUILD_DIR_MODULE)/private/
+$(BUILD_DIR_MODULE)/private/gauss_fermi.f90: ../private/gauss_fermi.f90 | $(BUILD_DIR_MODULE)/private/.
 	@cp $^ $@
 
-$(BUILD_DIR_MODULE)/private/helm.f90: ../private/helm.f90 | $(BUILD_DIR_MODULE)/private/
+$(BUILD_DIR_MODULE)/private/helm.f90: ../private/helm.f90 | $(BUILD_DIR_MODULE)/private/.
 	@cp $^ $@

eos/eosDT_builder/src/create_eos_files.f90

@@ -32,7 +32,7 @@ program create_eosDT_files
    real(dp) :: logQ_min = -10.09d0, logQ_max = 5.69d0, logRho_min = -14d0
    real(dp) :: logT_max
 
-   integer, parameter :: version_number = 51 ! update this to force rebuilding of caches
+   integer, parameter :: version_number = 52 ! update this to force rebuilding of caches
    ! update min_version in eosDT_load_tables to force rebuild of data files
 
    integer :: ix, io_unit, ios, info, irad
@@ -97,7 +97,8 @@ subroutine Make_EoS_Files(Z_in, X_in, include_radiation)
                        opal_only = .false., scvh_only = .false., search_for_SCVH = .true.
 
       !opal_only = .true.
-      scvh_only = .true.
+      opal_scvh_only = .true.
+      scvh_only = .false.
 
       dir = 'data' ! where to put the new data files
       io_unit = 40
@@ -178,7 +179,7 @@ subroutine Make_EoS_Files(Z_in, X_in, include_radiation)
                call do_stop('failed in helm_opal_scvh')
             end if
 
-            write (io_unit, '(f4.2,3(f10.5),7(1pe13.5),1(0pf9.5),4(0pf10.5),2(0pf11.5))') &
+            write (io_unit, '(f4.2,3(1x,f10.5),7(1x,1pe13.5),1x,0pf9.5,1x,f12.5,3(1x,1pe13.5),2(1x,0pf12.5))') &
                logT, logPgas, logE, logS, chiRho, chiT, Cp, Cv, dE_dRho, dS_dT, dS_dRho, mu, &
                log10(max(1d-99, free_e)), gamma1, gamma3, grad_ad, eta, HELM_fraction
          end do

eos/eosDT_builder/src/eos_regions_code.dek

@@ -1,4 +1,4 @@
-         if (logT >= logT1 .or. logT <= logT7 .or. logRho >= logRho1) then
+         if (logT >= logT1 .or. logT < logT7 .or. logRho >= logRho1) then
             iregion = pure_helm
          else if (logRho - 2*logT + 12 < logQmin) then
             iregion = pure_helm

eos/eosDT_builder/src/eos_regions_defs.dek

@@ -15,10 +15,11 @@
          double precision, parameter :: logT2_no_rad_default   = 7.7d0
                ! problems with HELM blend if logT2_no_rad_default is < logT1_no_rad_default
 
-         double precision, parameter :: logRho5 = -9.0d0
-         double precision, parameter :: logRho6 = -9.99d0
-         double precision, parameter :: logRho7 = -13.5d0
-         double precision, parameter :: logRho8 = -13.9d0
+         ! keep the generated table out of the warm low-density FreeEOS side
+         double precision, parameter :: logRho5 = -14.299d0
+         double precision, parameter :: logRho6 = -14.3d0
+         double precision, parameter :: logRho7 = -14.90d0
+         double precision, parameter :: logRho8 = -14.99d0
          double precision, parameter :: logT4   = 3.6d0
          double precision, parameter :: logT5   = 3.5d0
          ! 2.10 is lowest in SCVH; make mesa table use all of SCVH range at lowT end
@@ -33,7 +34,7 @@
          double precision, parameter :: logRho4 = logQ2 + 2*logT6 - 12
 
          double precision, parameter :: logT_HELM   = 3d0  ! smallest logT for HELM
-         double precision, parameter :: logQmin = -8.2d0
+         double precision, parameter :: logQmin = -8.05d0
          double precision, parameter :: logQmax = logQ1
 
-         double precision :: alfa, beta, c_dx, c_dy, logRho_lo, logRho_hi, logT1, logT2
+         double precision :: alfa, beta, c_dx, c_dy, logRho_lo, logRho_hi, logT1, logT2

eos/eosDT_builder/src/helm_opal_scvh_driver.f90

@@ -151,7 +151,7 @@ subroutine helm_opal_scvh(helm_only, opal_scvh_only, opal_only, scvh_only, &
       use helm
       use opal_scvh_driver
       use utils_lib, only: is_bad_num
-      use const_def, only: crad, ln10
+      use const_def, only: crad, kerg, amu
 
 !..logT and logRho are log10's of temp and den
       implicit none
@@ -371,14 +371,18 @@ subroutine helm_opal_scvh(helm_only, opal_scvh_only, opal_only, scvh_only, &
 
       subroutine check_results
          include 'formats'
+
+         ! Only reject non-finite or unphysical OPAL/SCVH support points.
+         ! Falling back to HELM for finite pressure-ionization structure
+         ! causes noisy mu features in regenerated tables.
          if (is_bad_num(logPgas_opalscvh) .or. is_bad_num(logE_opalscvh) .or. is_bad_num(logS_opalscvh) .or. &
              is_bad_num(chiRho_opalscvh) .or. is_bad_num(chiT_opalscvh) .or. is_bad_num(Cp_opalscvh) .or. &
              is_bad_num(Cv_opalscvh) .or. is_bad_num(dE_dRho_opalscvh) .or. is_bad_num(dS_dT_opalscvh) .or. &
              is_bad_num(dS_dRho_opalscvh) .or. is_bad_num(mu_opalscvh) .or. is_bad_num(logNe_opalscvh) .or. &
              is_bad_num(gamma1_opalscvh) .or. is_bad_num(gamma3_opalscvh) .or. is_bad_num(grad_ad_opalscvh) .or. &
              is_bad_num(eta_opalscvh) .or. Cv_opalscvh <= 0d0 .or. logPgas_opalscvh <= -50d0 .or. &
              logE_opalscvh <= -50d0 .or. logS_opalscvh <= -50d0 .or. logPgas_opalscvh <= -50d0 .or. &
-             gamma1_opalscvh <= 1d0 .or. gamma1_opalscvh >= 2d0 .or. grad_ad_opalscvh >= 50d0) then
+             gamma1_opalscvh <= 0d0) then
             ierr = -1
             !write(*,*) 'bail to HELM for logT logRho logQ', logT, logRho, logRho - 2*logT + 12
             if (.false.) then
@@ -414,19 +418,14 @@ subroutine get_helmeos(include_radiation)
          call helmeos2(temp, logT, den, logRho, abar, zbar, 1d6, 1d3, &
                        helm_res, clip_to_table_boundaries, .true., include_elec_pos, off_table, ierr)
          if (ierr /= 0) then
-            write (*, *) 'failed in helmeos2'
-            write (*, 1) 'temp', temp
-            write (*, 1) 'logT', logT
-            write (*, 1) 'den', den
-            write (*, 1) 'logRho', logRho
-            write (*, 1) 'Z', Z
-            write (*, 1) 'X', X
-            write (*, 1) 'abar', abar
-            write (*, 1) 'zbar', zbar
-            write (*, *) 'clip_to_table_boundaries', clip_to_table_boundaries
-            write (*, *) 'include_radiation', include_radiation
-            write (*, *) 'stop in get_helmeos'
-            stop 1
+            call helmeos2(temp, logT, den, logRho, abar, zbar, 1d6, 1d3, &
+                          helm_res, clip_to_table_boundaries, .true., .false., off_table, ierr)
+         end if
+         if (ierr /= 0) then
+            ! Last resort for rectangular eosDT support points outside HELM coverage.
+            call get_ideal_helmeos(include_radiation)
+            ierr = 0
+            return
          end if
 
          logPgas_helm = log10(helm_res(h_pgas))
@@ -460,6 +459,63 @@ subroutine get_helmeos(include_radiation)
 
       end subroutine get_helmeos
 
+      ! Used only when HELM cannot fill an eosDT support point.
+      subroutine get_ideal_helmeos(include_radiation)
+         logical, intent(in) :: include_radiation
+         double precision :: mu_ideal, Rgas, Pgas, Prad, P, egas, erad, ener, &
+            sgas, srad, entr, dpressdd, dpressdt, dedt, dedd, dsdt, dsdd, &
+            gamma3_minus1
+
+         mu_ideal = abar/(1d0 + zbar)
+         Rgas = kerg/(mu_ideal*amu)
+
+         Pgas = den*Rgas*temp
+         egas = 1.5d0*Rgas*temp
+         sgas = Rgas*(1.5d0*log(temp) - log(den) + 50d0)
+
+         Prad = 0d0
+         erad = 0d0
+         srad = 0d0
+         if (include_radiation) then
+            Prad = crad*temp**4/3d0
+            erad = 3d0*Prad/den
+            srad = 4d0*Prad/(den*temp)
+         end if
+
+         P = Pgas + Prad
+         ener = egas + erad
+         entr = sgas + srad
+
+         dpressdd = Rgas*temp
+         dpressdt = den*Rgas + 4d0*Prad/temp
+
+         dedt = 1.5d0*Rgas + 4d0*erad/temp
+         dedd = -erad/den
+
+         dsdt = 1.5d0*Rgas/temp + 3d0*srad/temp
+         dsdd = -(Rgas + srad)/den
+
+         logPgas_helm = log10(Pgas)
+         logE_helm = log10(ener)
+         logS_helm = log10(entr)
+
+         chiRho_helm = dpressdd*den/P
+         chiT_helm = dpressdt*temp/P
+         Cv_helm = dedt
+         Cp_helm = Cv_helm + P*chiT_helm**2/(den*temp*chiRho_helm)
+         gamma3_minus1 = dpressdt/(den*dedt)
+         gamma3_helm = 1d0 + gamma3_minus1
+         gamma1_helm = chiRho_helm + chiT_helm*gamma3_minus1
+         grad_ad_helm = gamma3_minus1/gamma1_helm
+         dE_dRho_helm = dedd
+         dS_dRho_helm = dsdd
+         dS_dT_helm = dsdt
+         logNe_helm = log10(max(1d-99, den*avo*zbar/abar))
+         mu_helm = mu_ideal
+         eta_helm = -20d0
+
+      end subroutine get_ideal_helmeos
+
    end subroutine helm_opal_scvh
 
 end module helm_opal_scvh_driver