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..75a5fb440 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -39,6 +39,26 @@ 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. + +At the low-density OPAL/SCVH edge, MESA now blends the OPAL/SCVH component to +the ideal EOS over :math:`-10 \leq \log_{10}\rho \leq -8` in the low-T +OPAL/SCVH region where MESA was previously falling back to HELM. + +Fixed the OPAL/SCVH blend between ``Z = 0.035`` and ``Z = 0.04``. This blend +was not working properly, so the EOS was not actually blending in Z over that +interval. For more details, see +:ref:`the known bugs entry ` and +`gh-997 `_. + +.. figure:: changelog_plots/eos_regions_x070_z002_compare.png + :alt: EOS regions in mesa-r26.4.1 and this release for X = 0.7, Z = 0.02 + + EOS regions for a solar-like composition in mesa-r26.4.1 and this release. + .. 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/changelog_plots/eos_regions_x070_z002_compare.png b/docs/source/changelog_plots/eos_regions_x070_z002_compare.png new file mode 100644 index 000000000..0b77d8b78 Binary files /dev/null and b/docs/source/changelog_plots/eos_regions_x070_z002_compare.png differ diff --git a/docs/source/eos/eos_regions_x000_z100.png b/docs/source/eos/eos_regions_x000_z100.png index 472b66223..715e6f898 100644 Binary files a/docs/source/eos/eos_regions_x000_z100.png and b/docs/source/eos/eos_regions_x000_z100.png differ diff --git a/docs/source/eos/eos_regions_x070_z002.png b/docs/source/eos/eos_regions_x070_z002.png index 58f0c7e04..6f6201ff2 100644 Binary files a/docs/source/eos/eos_regions_x070_z002.png and b/docs/source/eos/eos_regions_x070_z002.png differ diff --git a/docs/source/known_bugs.rst b/docs/source/known_bugs.rst index 64abefa1b..35332433b 100644 --- a/docs/source/known_bugs.rst +++ b/docs/source/known_bugs.rst @@ -13,6 +13,33 @@ issue, but it may not be complete. r26.4.1 ======= +.. _freeeos_fallback_log_units_bug: + +EOS: FreeEOS fallback rows and OPAL/SCVH blends +----------------------------------------------- + +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. + +This also exposed a low-density OPAL/SCVH edge where MESA was falling back to +HELM. The fix updates this region so OPAL/SCVH blends to the ideal EOS at low +density instead. + +Another OPAL/SCVH blend issue affected the high-Z transition between +``Z = 0.035`` and ``Z = 0.04``. In practice these regions were not being +blended in Z: OPAL/SCVH remained fully active below ``Z = 0.04`` wherever it +could be used, then switched discontinuously to the HELM/ideal layer at +``Z = 0.04``. This happened because the Z ramp multiplied the existing fallback +fraction, which was zero in pure OPAL/SCVH regions. The fix fades OPAL/SCVH +into the HELM/ideal layer across this Z interval. + +This is fixed after ``r26.4.1``. See +`gh-921 `_. Users who need the +patch directly can refer to `gh-997 `_. + .. _report_max_infall_inside_fe_core_bug: Controls: ``report_max_infall_inside_fe_core`` is ignored diff --git a/eos/defaults/eos.defaults b/eos/defaults/eos.defaults index 5373861c4..8eb5627b6 100644 --- a/eos/defaults/eos.defaults +++ b/eos/defaults/eos.defaults @@ -101,8 +101,8 @@ logRho_max_FreeEOS_hi = 10.1d0 logRho_max_FreeEOS_lo = 10.0d0 - logT_min_FreeEOS_hi = 3.1d0 - logT_min_FreeEOS_lo = 3.0d0 + logT_min_FreeEOS_hi = 3.12d0 + logT_min_FreeEOS_lo = 3.02d0 logT_max_FreeEOS_hi = 8.2d0 ! FreeEOS tabulations go up to here logT_max_FreeEOS_lo = 8.1d0 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..7f8d27703 100644 --- a/eos/eosDT_builder/src/scvh_core.f90 +++ b/eos/eosDT_builder/src/scvh_core.f90 @@ -402,6 +402,11 @@ 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 + info = -1 + 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 @@ -654,7 +659,6 @@ subroutine interpolate_scvh(include_radiation, search_for_SCVH, logT, logRho, T, si0t, si1t, si2t, si0mt, si1mt, si2mt, si0p, si1p, si2p, si0mp, si1mp, si2mp, logRho_test, & logP_guess, epslogP, epslogRho, x1, x3, y1, y3, dfdlogP, & z, fi_h(36), fi_he(36), w0t, w1t, w2t, w0mt, w1mt, w2mt, w0p, w1p, w2p, w0mp, w1mp, w2mp - include 'formats' info = 0 ipar => ipar_array @@ -936,10 +940,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))) diff --git a/eos/eosDT_data.tar.xz b/eos/eosDT_data.tar.xz index 2b27fba59..45f9b0c0a 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:c1c3062dbdfe4ca332ca12fe317c77d054e69906401523e68812d28dca690b16 +size 71063492 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/private/eosdt_eval.f90 b/eos/private/eosdt_eval.f90 index 607173a6f..9208bb7de 100644 --- a/eos/private/eosdt_eval.f90 +++ b/eos/private/eosdt_eval.f90 @@ -768,6 +768,12 @@ subroutine get_HELM_alfa( & ierr = 0 ht => eos_ht helm_blend_width = 0.1d0 + if (rq% use_FreeEOS .and. logT <= rq% logT_min_FreeEOS_lo) then + alfa = 1d0 + d_alfa_dlogRho = 0d0 + d_alfa_dlogT = 0d0 + return + end if bounds(1,1) = ht% logdlo bounds(1,2) = ht% logthi @@ -979,6 +985,7 @@ subroutine get_opal_scvh_for_eosdt( & rho, logRho, T, logT, remaining_fraction, & res, d_dlnd, d_dlnT, d_dxa, & skip, ierr) + use eosdt_support, only: Do_Blend integer, intent(in) :: handle logical, intent(in) :: dbg real(dp), intent(in) :: & @@ -992,26 +999,103 @@ subroutine get_opal_scvh_for_eosdt( & real(dp), intent(inout), dimension(nv, species) :: d_dxa logical, intent(out) :: skip integer, intent(out) :: ierr + type (EoS_General_Info), pointer :: rq + real(dp), parameter :: logRho_tail_ideal = -10d0 + real(dp), parameter :: logRho_tail_scvh = -8d0 + real(dp), dimension(nv) :: res_scvh, d_dlnd_scvh, d_dlnT_scvh, & + res_ideal, d_dlnd_ideal, d_dlnT_ideal + real(dp), dimension(nv, species) :: d_dxa_scvh, d_dxa_ideal + real(dp) :: tail_alfa, d_tail_alfa_dlogRho, d_tail_alfa_dlogT, & + logT_tail_off + logical :: skip_scvh, skip_ideal + + ierr = 0 + skip = .false. + rq => eos_handles(handle) + logT_tail_off = rq% logT_min_FreeEOS_hi + + if (logRho <= logRho_tail_ideal .and. logT <= logT_tail_off) then + call get_ideal_for_eosdt( & + handle, dbg, Z, X, abar, zbar, & + species, chem_id, net_iso, xa, & + rho, logRho, T, logT, remaining_fraction, & + res, d_dlnd, d_dlnT, d_dxa, & + skip, ierr) + return + end if + call get1_for_eosdt( & handle, eosdt_OPAL_SCVH, dbg, & Z, X, abar, zbar, & species, chem_id, net_iso, xa, & rho, logRho, T, logT, remaining_fraction, & - res, d_dlnd, d_dlnT, d_dxa, & - skip, ierr) + res_scvh, d_dlnd_scvh, d_dlnT_scvh, d_dxa_scvh, & + skip_scvh, ierr) + if (ierr /= 0) return + if (skip_scvh) then + skip = .true. + return + end if - ! zero phase information - res(i_phase:i_latent_ddlnRho) = 0d0 - d_dlnT(i_phase:i_latent_ddlnRho) = 0d0 - d_dlnd(i_phase:i_latent_ddlnRho) = 0d0 + call mark_opal_scvh_result(res_scvh, d_dlnd_scvh, d_dlnT_scvh) - ! zero all components - res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0 - d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0 - d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0 + if (logRho >= logRho_tail_scvh .or. logT >= logT_tail_off) then + res = res_scvh + d_dlnd = d_dlnd_scvh + d_dlnT = d_dlnT_scvh + d_dxa = d_dxa_scvh + skip = .false. + return + end if - ! mark this one - res(i_frac_OPAL_SCVH) = 1.0d0 + call get_ideal_for_eosdt( & + handle, dbg, Z, X, abar, zbar, & + species, chem_id, net_iso, xa, & + rho, logRho, T, logT, remaining_fraction, & + res_ideal, d_dlnd_ideal, d_dlnT_ideal, d_dxa_ideal, & + skip_ideal, ierr) + if (ierr /= 0) return + if (skip_ideal) then + skip = .true. + return + end if + + tail_alfa = (logRho_tail_scvh - logRho)/(logRho_tail_scvh - logRho_tail_ideal) + if (tail_alfa <= 0d0) then + tail_alfa = 0d0 + d_tail_alfa_dlogRho = 0d0 + else if (tail_alfa >= 1d0) then + tail_alfa = 1d0 + d_tail_alfa_dlogRho = 0d0 + else + d_tail_alfa_dlogRho = -1d0/(logRho_tail_scvh - logRho_tail_ideal) + end if + d_tail_alfa_dlogT = 0d0 + + call Do_Blend( & + rq, species, rho, logRho, T, logT, & + tail_alfa, d_tail_alfa_dlogT, d_tail_alfa_dlogRho, .false., & + res_ideal, d_dlnd_ideal, d_dlnT_ideal, d_dxa_ideal, & + res_scvh, d_dlnd_scvh, d_dlnT_scvh, d_dxa_scvh, & + res, d_dlnd, d_dlnT, d_dxa) + skip = .false. + + contains + + subroutine mark_opal_scvh_result(res_in, d_dlnd_in, d_dlnT_in) + real(dp), intent(inout), dimension(nv) :: res_in, d_dlnd_in, d_dlnT_in + + res_in(i_phase:i_latent_ddlnRho) = 0d0 + d_dlnT_in(i_phase:i_latent_ddlnRho) = 0d0 + d_dlnd_in(i_phase:i_latent_ddlnRho) = 0d0 + + res_in(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0 + d_dlnd_in(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0 + d_dlnT_in(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0 + + res_in(i_frac_OPAL_SCVH) = 1.0d0 + + end subroutine mark_opal_scvh_result end subroutine get_opal_scvh_for_eosdt @@ -1448,11 +1532,11 @@ subroutine set_alfa_and_partials ! alfa = fraction other end if if (Z > Z_no_HELM .and. Z < Z_all_HELM .and. alfa < 1d0) then - ! reduce alfa to reduce the HELM fraction + ! Increase alfa with Z to fade OPAL/SCVH into the next EOS. zfactor = (Z - Z_no_HELM)/(Z_all_HELM - Z_no_HELM) - alfa = alfa*zfactor - d_alfa_dlogRho = d_alfa_dlogRho*zfactor - d_alfa_dlogT = d_alfa_dlogT*zfactor + alfa = zfactor + (1d0 - zfactor)*alfa + d_alfa_dlogRho = (1d0 - zfactor)*d_alfa_dlogRho + d_alfa_dlogT = (1d0 - zfactor)*d_alfa_dlogT end if end subroutine set_alfa_and_partials 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