diff --git a/atm/test/test_output b/atm/test/test_output index 8718566aa..3560d0d46 100644 --- a/atm/test/test_output +++ b/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 diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 34fc942a7..c6261f9f6 100644 --- a/docs/source/changelog.rst +++ b/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``). diff --git a/docs/source/known_bugs.rst b/docs/source/known_bugs.rst index 64abefa1b..1ef5ecbcd 100644 --- a/docs/source/known_bugs.rst +++ b/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 diff --git a/eos/eosDT_builder/Makefile b/eos/eosDT_builder/Makefile index ad5637318..4cf05518d 100644 --- a/eos/eosDT_builder/Makefile +++ b/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 $^ $@ diff --git a/eos/eosDT_builder/src/create_eos_files.f90 b/eos/eosDT_builder/src/create_eos_files.f90 index 9a70f6e51..3d2064b83 100644 --- a/eos/eosDT_builder/src/create_eos_files.f90 +++ b/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 diff --git a/eos/eosDT_builder/src/eos_regions_code.dek b/eos/eosDT_builder/src/eos_regions_code.dek index d4d2b510e..325ea01b3 100644 --- a/eos/eosDT_builder/src/eos_regions_code.dek +++ b/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 diff --git a/eos/eosDT_builder/src/eos_regions_defs.dek b/eos/eosDT_builder/src/eos_regions_defs.dek index 2fe378634..2143f4c7d 100644 --- a/eos/eosDT_builder/src/eos_regions_defs.dek +++ b/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 \ No newline at end of file + double precision :: alfa, beta, c_dx, c_dy, logRho_lo, logRho_hi, logT1, logT2 diff --git a/eos/eosDT_builder/src/helm_opal_scvh_driver.f90 b/eos/eosDT_builder/src/helm_opal_scvh_driver.f90 index dfc210f51..a5cc26e2d 100644 --- a/eos/eosDT_builder/src/helm_opal_scvh_driver.f90 +++ b/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,6 +371,10 @@ 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. & @@ -378,7 +382,7 @@ subroutine check_results 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 diff --git a/eos/eosDT_builder/src/scvh_core.f90 b/eos/eosDT_builder/src/scvh_core.f90 index 831f38fdc..8586186fd 100644 --- a/eos/eosDT_builder/src/scvh_core.f90 +++ b/eos/eosDT_builder/src/scvh_core.f90 @@ -32,6 +32,8 @@ module scvh_core integer, parameter :: NlogPs = 99, NlogTs = 63 double precision, parameter :: dlogP = 0.20d0 double precision, parameter :: logP_min = -0.60d0 + ! Analytic continuation below the processed SCVH edge. + double precision, parameter :: logP_tail_min = -6.00d0 double precision, parameter :: logP_max = 19d0 double precision, parameter :: dlogT = 0.08d0 double precision, parameter :: logT_min = 2.10d0 @@ -51,6 +53,8 @@ module scvh_core type(scvh_info), target ::compv1_h_si, compv2_h_si, denlog_h_si, slog_h_si, ulog_h_si, & compv1_he_si, compv2_he_si, denlog_he_si, slog_he_si, ulog_he_si + double precision :: slog_edge_h(NlogTs), slog_edge_he(NlogTs) + contains subroutine setup_scvh(data_dir) @@ -122,6 +126,10 @@ subroutine read_data_for_scvh(data_dir) call read_file_for_scvh(data_dir, 'scvh/he_tab.asc.data', & num_he_pts, NlogPs, NlogTs, compv1_he_si, compv2_he_si, denlog_he_si, slog_he_si, ulog_he_si) + ! Save the restored entropy edge before padding the rectangular spline table. + slog_edge_h(:) = slog_h_si%Z(1, :) + slog_edge_he(:) = slog_he_si%Z(1, :) + if (.false.) then call write_file_for_scvh(data_dir, 'scvh/h_tab.new.data', & num_h_pts, NlogPs, NlogTs, compv1_h_si, compv2_h_si, denlog_h_si, slog_h_si, ulog_h_si) @@ -402,6 +410,12 @@ subroutine interp_vals_bicub(only_densities, search_for_SCVH, den_h, ener_h, ent max_logP = logPs(kP) end if + if (logP < logP_min) then + ! Below the processed edge, use an edge-matched ideal-gas tail. + call interp_vals_lowP_tail(only_densities) + return + end if + den_h = 10.0**interp_value(num_h_pts, denlog_h_si); if (info /= 0) return dddpress_ct_h = d_dlogP dddt_cp_h = d_dlogT @@ -518,6 +532,184 @@ subroutine interp_vals_bicub(only_densities, search_for_SCVH, den_h, ener_h, ent contains + subroutine interp_vals_lowP_tail(only_densities_arg) + logical, intent(in) :: only_densities_arg + + call tail_species(num_h_pts, compv1_h_si, compv2_h_si, denlog_h_si, slog_h_si, ulog_h_si, slog_edge_h, & + den_h, ener_h, entr_h, dddt_cp_h, dddpress_ct_h, dsdt_cp_h, dsdpress_ct_h, & + xnh2, dxnh2_dlogT, dxnh2_dlogP, xnh, dxnh_dlogT, dxnh_dlogP) + if (info /= 0) return + + call tail_species(num_he_pts, compv1_he_si, compv2_he_si, denlog_he_si, slog_he_si, ulog_he_si, slog_edge_he, & + den_he, ener_he, entr_he, dddt_cp_he, dddpress_ct_he, dsdt_cp_he, dsdpress_ct_he, & + xnhe, dxnhe_dlogT, dxnhe_dlogP, xnhep, dxnhep_dlogT, dxnhep_dlogP) + if (info /= 0) return + + if (only_densities_arg) return + + call limit_pair(xnh2, dxnh2_dlogP, dxnh2_dlogT, xnh, dxnh_dlogP, dxnh_dlogT) + call limit_pair(xnhe, dxnhe_dlogP, dxnhe_dlogT, xnhep, dxnhep_dlogP, dxnhep_dlogT) + + end subroutine interp_vals_lowP_tail + + subroutine tail_species(n, comp1_si, comp2_si, denlog_si, slog_si, ulog_si, slog_edge, & + den, ener, entr, dddt_cp, dddpress_ct, dsdt_cp, dsdpress_ct, & + x1, dx1_dlogT, dx1_dlogP, x2, dx2_dlogT, dx2_dlogP) + integer, intent(in) :: n + type(scvh_info) :: comp1_si, comp2_si, denlog_si, slog_si, ulog_si + double precision, intent(in) :: slog_edge(NlogTs) + double precision, intent(out) :: den, ener, entr, dddt_cp, dddpress_ct, dsdt_cp, dsdpress_ct, & + x1, dx1_dlogT, dx1_dlogP, x2, dx2_dlogT, dx2_dlogP + + double precision :: logden0, dlogden0_dlogP, dlogden0_dlogT, logden, dlogden_dlogP, dlogden_dlogT, & + logs0, dlogs0_dlogP, dlogs0_dlogT, logs, dlogs_dlogP, dlogs_dlogT, & + logu0, dlogu0_dlogP, dlogu0_dlogT, logu, dlogu_dlogP, dlogu_dlogT, & + x10, dx10_dlogP, dx10_dlogT, x20, dx20_dlogP, dx20_dlogT, & + logP_decay, S0, S, Rgas, dRgas_dlogT, dS_dlogT, ideal_logs, ideal_dlogs_dlogP, ideal_dlogs_dlogT + + ! Anchor the tail to the restored processed table at logP_min. + call edge_value(n, denlog_si, logden0, dlogden0_dlogP, dlogden0_dlogT) + if (info /= 0) return + + logP_decay = max(logP_tail_min, logP_min - dlogP) + + call edge_value(n, comp1_si, x10, dx10_dlogP, dx10_dlogT) + if (info /= 0) return + call edge_value(n, comp2_si, x20, dx20_dlogP, dx20_dlogT) + if (info /= 0) return + call edge_value(n, ulog_si, logu0, dlogu0_dlogP, dlogu0_dlogT) + if (info /= 0) return + call edge_value(n, slog_si, logs0, dlogs0_dlogP, dlogs0_dlogT) + if (info /= 0) return + dlogs0_dlogT = edge_dlogT(slog_si, slog_edge) + + call tail_value(dx10_dlogP, x10, 0d0, dx10_dlogT, 0d0, logP_decay, x1, dx1_dlogP, dx1_dlogT) + call tail_value(dx20_dlogP, x20, 0d0, dx20_dlogT, 0d0, logP_decay, x2, dx2_dlogP, dx2_dlogT) + + call tail_value(dlogden0_dlogP, logden0 + logP - logP_min, & + 1d0, -1d0, 1d0, logP_decay, logden, dlogden_dlogP, dlogden_dlogT) + + call tail_value(dlogu0_dlogP, logu0, 0d0, dlogu0_dlogT, 0d0, logP_decay, & + logu, dlogu_dlogP, dlogu_dlogT) + + S0 = 10d0**logs0 + Rgas = 10d0**(logP_min - logden0 - logT) + S = S0 + Rgas*log(10d0)*(logP_min - logP) + dRgas_dlogT = 0d0 + dS_dlogT = log(10d0)*S0*dlogs0_dlogT + dRgas_dlogT*log(10d0)*(logP_min - logP) + ideal_logs = log10(S) + ideal_dlogs_dlogP = -Rgas/S + ideal_dlogs_dlogT = dS_dlogT/(log(10d0)*S) + + call tail_value(dlogs0_dlogP, ideal_logs, ideal_dlogs_dlogP, ideal_dlogs_dlogT, & + -Rgas/S0, logP_decay, logs, dlogs_dlogP, dlogs_dlogT) + + den = 10d0**logden + dddpress_ct = dlogden_dlogP + dddt_cp = dlogden_dlogT + + ener = 10d0**logu + entr = 10d0**logs + dsdpress_ct = dlogs_dlogP + dsdt_cp = dlogs_dlogT + + end subroutine tail_species + + subroutine edge_value(n, si, val, dval_dlogP, dval_dlogT) + integer, intent(in) :: n + type(scvh_info) :: si + double precision, intent(out) :: val, dval_dlogP, dval_dlogT + + XI(1) = logP_min + val = interp_value(n, si) + if (info /= 0) return + dval_dlogP = d_dlogP + dval_dlogT = d_dlogT + + end subroutine edge_value + + double precision function edge_dlogT(si, edge) + type(scvh_info) :: si + double precision, intent(in) :: edge(NlogTs) + integer :: j + double precision :: v0, v1 + + j = locate_logT(logT) + if (j < 1) j = 1 + if (j >= NlogTs) j = NlogTs - 1 + + v0 = edge(j) + v1 = edge(j + 1) + if (v0 == 0d0 .or. v1 == 0d0) then + v0 = si%Z(1, j) + v1 = si%Z(1, j + 1) + end if + edge_dlogT = (v1 - v0)/(logTs(j + 1) - logTs(j)) + + end function edge_dlogT + + subroutine tail_value(edge_dlogP, ideal_val, ideal_dlogP, ideal_dlogT, & + ideal_edge_dlogP, logP_decay, val, dval_dlogP, dval_dlogT) + double precision, intent(in) :: edge_dlogP, ideal_val, ideal_dlogP, ideal_dlogT, ideal_edge_dlogP, logP_decay + double precision, intent(out) :: val, dval_dlogP, dval_dlogT + double precision :: L, x, g, dg_dx, M + + if (logP <= logP_decay) then + val = ideal_val + dval_dlogP = ideal_dlogP + dval_dlogT = ideal_dlogT + return + end if + + L = logP_min - logP_decay + x = (logP - logP_decay)/L + ! C1 correction: zero value at both ends, zero slope at the + ! ideal end, and unit slope at the SCVH edge. + g = x**3 - x**2 + dg_dx = 3d0*x**2 - 2d0*x + M = (edge_dlogP - ideal_edge_dlogP)*L + + val = ideal_val + M*g + dval_dlogP = ideal_dlogP + M*dg_dx/L + dval_dlogT = ideal_dlogT + + end subroutine tail_value + + subroutine limit_pair(x1, dx1_dlogP, dx1_dlogT, x2, dx2_dlogP, dx2_dlogT) + double precision, intent(inout) :: x1, dx1_dlogP, dx1_dlogT, x2, dx2_dlogP, dx2_dlogT + + call limit_fraction(x1, dx1_dlogP, dx1_dlogT) + call limit_fraction(x2, dx2_dlogP, dx2_dlogT) + + if (x1 + x2 > 1d0) then + if (x1 > x2) then + x1 = 1d0 - x2 + dx1_dlogP = -dx2_dlogP + dx1_dlogT = -dx2_dlogT + else + x2 = 1d0 - x1 + dx2_dlogP = -dx1_dlogP + dx2_dlogT = -dx1_dlogT + end if + end if + + end subroutine limit_pair + + subroutine limit_fraction(x, dx_dlogP, dx_dlogT) + double precision, intent(inout) :: x, dx_dlogP, dx_dlogT + + if (x > 1d0) then + x = 1d0 + dx_dlogP = 0d0 + dx_dlogT = 0d0 + else if (x < 0d0) then + x = 0d0 + dx_dlogP = 0d0 + dx_dlogT = 0d0 + end if + + end subroutine limit_fraction + double precision function interp_value(n, si) use interp_2d_lib_db, only: interp_mkbicub_db, interp_evbicub_db use num_lib @@ -609,7 +801,7 @@ subroutine interpolate_scvh(include_radiation, search_for_SCVH, logT, logRho, T, double precision :: logP, inv_Rho, inv_T, inv_P, small_value, Rho_min parameter(small_value=1.0d-16) - parameter(Rho_min=1.0d-10) + parameter(Rho_min=1.0d-15) !..for hydrogen double precision :: den_h, ener_h, entr_h, dddt_cp_h, dddpress_ct_h, & @@ -936,10 +1128,10 @@ subroutine interpolate_scvh(include_radiation, search_for_SCVH, logT, logRho, T, chiT = dpressdt*T/P gamma3 = 1 + dpressdt/(Rho*dedt) - grad_ad = dtdpress_cs_hhe/(T*inv_P) - gamma1 = (gamma3 - 1)/grad_ad ! C&G 9.88 & 9.89 - Cp = Cv + P*chiT**2/(Rho*T*chiRho) ! C&G 9.86 + ! Use total derivatives here since radiation may have been added above. + gamma1 = chiRho + chiT*(gamma3 - 1) ! C&G 9.88 + grad_ad = (gamma3 - 1)/gamma1 ! C&G 9.89 xnhp = max(0d0, min(1d0, 1 - (xnh2 + xnh))) xnhepp = max(0d0, min(1d0, 1 - (xnhep + xnhe))) @@ -1024,7 +1216,7 @@ double precision function fscvh(logP, dfdlogP, lrpar, rpar, lipar, ipar, ierr) dxnhe_dlogP, xnhep, dxnhep_dlogT, dxnhep_dlogP, xmassh1, xmasshe4 logical :: only_densities, search_for_SCVH include 'formats' - if (logP < logP_min .or. logP > logP_max) then + if (logP < logP_tail_min .or. logP > logP_max) then ierr = -1 if (dbg) write (*, *) if (dbg) write (*, 1) 'logP guess out of bounds', logP diff --git a/eos/eosDT_data.tar.xz b/eos/eosDT_data.tar.xz index 2b27fba59..0dc6ee602 100644 --- a/eos/eosDT_data.tar.xz +++ b/eos/eosDT_data.tar.xz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3992b7e130d01252a49cb4269d1de379b5fe289d6781839d087164c32cb2f674 -size 160326344 +oid sha256:de1a9e4c2267c0a782477a026b1879cc1c1e01f0745e422d7dde7997301c8215 +size 73954276 diff --git a/eos/eosFreeEOS_builder/inlist_co b/eos/eosFreeEOS_builder/inlist_co index 51ef9c8e6..8c1e8e58b 100644 --- a/eos/eosFreeEOS_builder/inlist_co +++ b/eos/eosFreeEOS_builder/inlist_co @@ -4,7 +4,7 @@ mass_list = 'mass_frac_CO.txt' !MESA/eos table version number - table_version = 51 + table_version = 52 !FreeEOS option suite, see FreeEOS README for details eos_version = 1 diff --git a/eos/eosFreeEOS_builder/inlist_solar b/eos/eosFreeEOS_builder/inlist_solar index e95b9089d..ac47f1e11 100644 --- a/eos/eosFreeEOS_builder/inlist_solar +++ b/eos/eosFreeEOS_builder/inlist_solar @@ -4,7 +4,7 @@ mass_list = 'mass_frac_solar.txt' !MESA/eos table version number - table_version = 51 + table_version = 52 !FreeEOS option suite, see FreeEOS README for details eos_version = 1 diff --git a/eos/eosFreeEOS_builder/src/free_eos_table.f90 b/eos/eosFreeEOS_builder/src/free_eos_table.f90 index bba4d9125..9742c53c8 100644 --- a/eos/eosFreeEOS_builder/src/free_eos_table.f90 +++ b/eos/eosFreeEOS_builder/src/free_eos_table.f90 @@ -248,7 +248,7 @@ subroutine read_namelist ! set defaults mass_list = 'mass_frac.txt' - table_version = 51 + table_version = 52 eos_version = 1 !set T range @@ -564,7 +564,12 @@ subroutine mesa_eos_eval( logRho0, logT, mass_Frac, eos_result) end if eos_result(1:num_eos_basic_results) = res - eos_result(i_lnRho) = logRho + + ! FreeEOS data files store these logarithms in base 10. + eos_result(i_lnPgas) = res(i_lnPgas)/ln10 + eos_result(i_lnE) = res(i_lnE)/ln10 + eos_result(i_lnS) = res(i_lnS)/ln10 + eos_result(i_lnRho) = log10Rho eos_result(i_dpe) = 0._dp eos_result(i_dsp) = 0._dp eos_result(i_dse) = 0._dp diff --git a/eos/eosFreeEOS_data.tar.xz b/eos/eosFreeEOS_data.tar.xz index 7d93c28d9..00fe48a88 100644 --- a/eos/eosFreeEOS_data.tar.xz +++ b/eos/eosFreeEOS_data.tar.xz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:fb7aed3e88b3ed295630f3c3ec8133c19fd3b0cce0580011474f1e5510ee47c1 -size 1386142312 +oid sha256:d927560ce390f2f7f534fbc154ec5122bf37fff01b2d3d2668ba29509d140bec +size 1379025232 diff --git a/eos/plotter/src/eos_plotter.f90 b/eos/plotter/src/eos_plotter.f90 index e1bb97825..94c22e807 100644 --- a/eos/plotter/src/eos_plotter.f90 +++ b/eos/plotter/src/eos_plotter.f90 @@ -4,6 +4,7 @@ program eos_plotter use eos_lib, only: eos_ptr, eosDT_get use chem_def use chem_lib + use const_def, only: crad, ln10 use const_lib, only: const_init use math_lib use num_lib, only : dfridr @@ -402,8 +403,8 @@ program eos_plotter end if - if (only_blend_regions .and. .not. in_eos_blend(res)) then - call set_nan(res1) + if (only_blend_regions) then + if (.not. in_eos_blend(res)) call set_nan(res1) end if write(iounit,*) kval, jval, res1 diff --git a/eos/test/test_output b/eos/test/test_output index 0ea2ff6a8..a0db323bf 100644 --- a/eos/test/test_output +++ b/eos/test/test_output @@ -36,15 +36,15 @@ eosDT_get P 2.3440634573327629E+10 - E 3.2446133432181909E+11 + E 3.2446106398564081E+11 S 2.8977788032200122E+08 - mu 4.0000007737813874E+00 + mu 4.0000007738172707E+00 lnfree_e -9.2103403719761818E+00 eta -2.0000000000000000E+01 - grad_ad 3.9197054769106937E-01 + grad_ad 3.9089425291681118E-01 chiRho 1.1192461048209497E+00 - chiT 9.4195159087688318E-01 - gamma1 1.7666744682116600E+00 + chiT 9.4195158827893422E-01 + gamma1 1.7715322485086398E+00 gamma3 1.6924793542208876E+00 eta -2.0000000000000000E+01 @@ -60,15 +60,15 @@ eosDT_get P 3.2383932609950977E+11 - E 9.4473235426827480E+12 + E 9.4473235426977852E+12 S 4.7301040764836609E+08 mu 1.9473573566866675E+00 - lnfree_e -1.3207182576098033E+00 + lnfree_e -1.3207366848530466E+00 eta -2.2140136870353762E+00 - grad_ad 2.5815798177346022E-01 - chiRho 8.8670604922952434E-01 + grad_ad 2.5815521408865183E-01 + chiRho 8.8670604663157537E-01 chiT 1.4448369312027372E+00 - gamma1 1.4141864668291979E+00 + gamma1 1.4141864011682150E+00 gamma3 1.3650780037467394E+00 eta -2.2140136870353762E+00 @@ -134,17 +134,17 @@ eosDT_get P 6.9084070757562546E+10 - E 1.8861617688433232E+12 + E 1.8861616560131572E+12 S 9.1014778380461836E+08 mu 1.3477398922872359E+00 - lnfree_e -7.7424749221373501E+00 - eta -1.6593753447884108E+01 - grad_ad 2.2741113591255283E-01 + lnfree_e -7.7425182669077861E+00 + eta -1.6593956642527232E+01 + grad_ad 2.2668862192983433E-01 chiRho 1.2436975923257054E+00 chiT 1.1782005930283301E+00 - gamma1 1.6915540427454960E+00 + gamma1 1.6969231782299208E+00 gamma3 1.3846744062743606E+00 - eta -1.6593753447884108E+01 + eta -1.6593956642527232E+01 X 1.0000000000000000E+00 @@ -158,15 +158,15 @@ eosDT_get P 1.2385332673497964E+12 - E 2.4388385774716191E+13 + E 2.4388385774571582E+13 S 1.4600490057583487E+09 mu 7.2181713563159555E-01 - lnfree_e -3.7774616224178065E-01 + lnfree_e -3.7774615138682827E-01 eta -1.2132438790845224E+00 - grad_ad 3.6117444052541564E-01 + grad_ad 3.6117916013083462E-01 chiRho 9.6642435620801737E-01 chiT 1.2770934029069718E+00 - gamma1 1.7938621097381113E+00 + gamma1 1.7938621097381100E+00 gamma3 1.6479051986892039E+00 eta -1.2132438790845224E+00 @@ -233,16 +233,16 @@ eosDT_get P 5.5654385291979820E+10 E 1.4844253665903733E+12 - S 7.6542994034342134E+08 + S 7.6542994035226464E+08 mu 2.0883727240864358E+00 - lnfree_e -8.2387144895465347E+00 - eta -1.8507079502615483E+01 - grad_ad 2.3708074718722938E-01 + lnfree_e -8.2387187722941917E+00 + eta -1.8507086026913143E+01 + grad_ad 2.3840400741225612E-01 chiRho 1.1951483607913123E+00 chiT 1.1869816367434431E+00 - gamma1 1.6762994130699729E+00 - gamma3 1.3973918063946891E+00 - eta -1.8507079502615483E+01 + gamma1 1.6669497128926525E+00 + gamma3 1.3973918058499728E+00 + eta -1.8507086026913143E+01 X 7.2750000000000004E-01 @@ -259,14 +259,14 @@ E 2.0332958756982102E+13 S 1.2178932952350833E+09 mu 9.2878901308633832E-01 - lnfree_e -5.3634909173875978E-01 - eta -1.3856884888818266E+00 - grad_ad 3.5162128251300456E-01 - chiRho 9.4786480115446770E-01 + lnfree_e -5.3634919068986409E-01 + eta -1.3856884879642708E+00 + grad_ad 3.5162275397339221E-01 + chiRho 9.4786480310326504E-01 chiT 1.2892278740569716E+00 - gamma1 1.7337719819494843E+00 + gamma1 1.7337719819487181E+00 gamma3 1.6095200113685473E+00 - eta -1.3856884888818266E+00 + eta -1.3856884879642708E+00 X 7.2750000000000004E-01 @@ -331,16 +331,16 @@ eosDT_get P 5.4349014621238190E+10 E 1.4432353712570466E+12 - S 7.4976378082837462E+08 + S 7.4976378083459973E+08 mu 2.1626082904518826E+00 - lnfree_e -8.2922278715096169E+00 - eta -1.8668741913261101E+01 - grad_ad 2.3836728669019633E-01 + lnfree_e -8.2922336208391982E+00 + eta -1.8668751839028083E+01 + grad_ad 2.3991422064089690E-01 chiRho 1.1904706709604147E+00 chiT 1.1866072226781790E+00 - gamma1 1.6752356524575822E+00 - gamma3 1.3992923413116880E+00 - eta -1.8668741913261101E+01 + gamma1 1.6643829175293345E+00 + gamma3 1.3992923407224072E+00 + eta -1.8668751839028083E+01 X 6.9999999999999996E-01 @@ -357,14 +357,14 @@ E 1.9942428832396027E+13 S 1.1939413753703966E+09 mu 9.5100508701167052E-01 - lnfree_e -5.5338792437867279E-01 - eta -1.4041113609979621E+00 - grad_ad 3.5011286661216673E-01 + lnfree_e -5.5338813597046310E-01 + eta -1.4041113609972320E+00 + grad_ad 3.5011726572624136E-01 chiRho 9.4570552362766236E-01 chiT 1.2902790434790619E+00 - gamma1 1.7248612299611938E+00 + gamma1 1.7248612299592598E+00 gamma3 1.6037880253708852E+00 - eta -1.4041113609979621E+00 + eta -1.4041113609972320E+00 X 6.9999999999999996E-01 @@ -431,26 +431,26 @@ logT 3.0000000000000000D+00 logPgas 4.2000000000000002D+00 - dlnRho_dlnPgas_c_T 9.9998987532966155D-01 - dlnRho_dlnT_c_Pgas -9.9975769898338618D-01 + dlnRho_dlnPgas_c_T 9.9998987532991812D-01 + dlnRho_dlnT_c_Pgas -9.9975769898364186D-01 P 1.5848934446522200D+04 - E 8.1826907094327927D+10 - S 8.9700174679926920D+08 + E 8.1826908401699203D+10 + S 8.9700174541741323D+08 mu 1.9000000000000001D+00 lnfree_e -9.2103403719761818D+00 eta -2.0000000000000000D+01 - grad_ad 2.9309721606956174D-01 - chiRho 1.0000099656494308D+00 - chiT 9.9976829870535666D-01 - Cp 1.2008175123022686D+08 - Cv 8.4993532780778155D+07 - dE_dRho 1.7987628811511891D+13 - dS_dT 8.4993532780778158D+04 - dS_dRho -7.7756530196169234D+13 - gamma1 1.4078416959964035D+00 - gamma3 1.4125438954496472D+00 + grad_ad 2.9212307792273379D-01 + chiRho 1.0000099656491741D+00 + chiT 9.9976829870535588D-01 + Cp 1.2008175049396412D+08 + Cv 8.4993566346818984D+07 + dE_dRho 1.7986469344859305D+13 + dS_dT 8.4993566346818989D+04 + dS_dRho -7.7756527275522500D+13 + gamma1 1.4124591784049583D+00 + gamma3 1.4125438954496474D+00 phase 0.0000000000000000D+00 latent_ddlnT 0.0000000000000000D+00 latent_ddlnRho 0.0000000000000000D+00