Merge pull request #615 from MESAHub/remove_HELM_floor

Remove HELM floor, fall back to ideal instead

Author: evbauer

docs/source/changelog.rst

@@ -24,7 +24,7 @@ pressure with respect to other variables.
 run_star_extras
 ~~~~~~~~~~~~~~~
 
-Previously we had logic to determine if a extra history value should be saved
+Previously we had logic to determine if an extra history value should be saved
 as an int or a float (users can only provide data as a float). This was error
 prone, so now we save extra history values as floats.
 
@@ -96,7 +96,7 @@ The phase diagram for O/Ne separation comes from
 `Blouin & Daligault (2021b) <https://ui.adsabs.harvard.edu/abs/2021ApJ...919...87B/abstract>`_.
 
 Massive Star test_suite Updates
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The ``12Msun_pre_ms_to_core_collapse``, ``20Msun_pre_ms_to_core_collapse``, and ``zams_to_cc_80``
 test_suites have each been updated and now fully evolve models at solar metallicities to core
@@ -123,6 +123,12 @@ The ``fixed_Teff``, ``fixed_Tsurf``, ``fixed_Psurf``,  and ``fixed_Psurf_and_Tsu
 atmosphere options were removed in r15140. We have reimplemented them although we
 caution users that their implementation could conflict with ``mlt_option = 'TDC'``.
 
+The EOS coverage regions have been updated to fall back to ideal gas in a region
+previously covered by HELM where it returned unphysical floor values of ``1e-20``
+for pressure, energy, and entropy. The most up-to-date EOS coverage plots can
+be found in the EOS documentation: :ref:`eos/overview:Overview of eos module`.
+
+
 Changes in r23.05.1
 ===================
 

docs/source/eos/eos_regions_x000_z100.png

@@ -43,7 +43,6 @@ def parse(fname):
 cbar_max = None
 
 pcol = ax.pcolormesh(Xran, Yran, eosDT[...,2], shading='nearest', cmap=cmap, vmin=cbar_min, vmax=cbar_max)
-pcol.set_edgecolor('face')
 cax = fig.colorbar(pcol, extend='both')
 cax.set_label('')
 

docs/source/eos/eos_regions_x070_z002.png

@@ -26,14 +26,22 @@ def parse(fname):
 # overwrite with fancier labels
 xlabel = r'$\log(\rho/{\rm g\,cm^{-3}})$'
 ylabel = r'$\log(T/{\rm K})$'
-title = r'MESA EOS Regions ($X=0.7$, $Z=0.02$)'
 
 eosDT, Yran, Xran = parse('eos_plotter.dat')
 
 apjcolwidth=3.38
 # set up plot and labels
 #fig, ax = plt.subplots(figsize=(apjcolwidth,apjcolwidth*4./5.)) # for paper figures
 fig, ax = plt.subplots(figsize=(5,4)) # for website pngs
+
+## for solar composition plot
+title = r'MESA EOS Regions ($X=0.7$, $Z=0.02$)'
+ax.text(4.0,2.3,'no coverage',fontsize='small')
+
+## for Z = 1 plot
+#title = r'MESA EOS Regions ($X=0.0$, $Z=1.0$)'
+#ax.text(4.8,2.3,'no coverage',fontsize='small')
+
 ax.set_title(title)
 ax.set_xlabel(xlabel)
 ax.set_ylabel(ylabel)
@@ -48,9 +56,9 @@ def parse(fname):
 cmap = colors.ListedColormap(my_colors)
 bounds = [-0.5, 0.5, 1.5, 2.5, 3.5, 5.5, 7.5]
 norm = colors.BoundaryNorm(bounds, cmap.N)
+cmap.set_under('white')
 
 pcol = ax.pcolormesh(Xran, Yran, eosDT[...,2], shading='nearest', cmap=cmap, norm=norm, rasterized=True)
-pcol.set_edgecolor('face')
 cax = fig.colorbar(pcol, ticks=[0, 1, 2, 3, 4.5, 6.5])
 cax.set_label('')
 cax.ax.tick_params(labelsize='x-small')

eos/plotter/plotter.py

@@ -257,6 +257,13 @@ subroutine helmeos2aux( &
          return
        end if
 
+       ! error out for very low pgas,egas,sgas (previously just set hard floor at this value)
+       if(pgas < 1d-20 .or. egas < 1d-20 .or. sgas < 1d-20) then
+          ierr = 1
+          if (dbg) write(*,*) 'failed in helm, sums too small'
+          return
+       end if
+
 !..compute the derivative quantities (cv, gamma1 ...etc)
        include 'helm_gammas.dek'
        if (ierr /= 0) then

eos/plotter/regions.py

@@ -1,4 +1,4 @@
-       pgas    = max(pion + pele + ppos + pcoul, 1d-20)
+       pgas    = pion + pele + ppos + pcoul
 
        dpgasdd  = dpiondd  + dpepdd  + dpcouldd 
        dpgasdt  = dpiondt  + dpepdt  + dpcouldt
@@ -15,7 +15,7 @@
        dpgasdaz = dpiondaz + dpepdaz + dpcouldaz 
        dpgasdzz = dpiondzz + dpepdzz + dpcouldzz 
 
-       egas    = max(eion + eele + epos + ecoul, 1d-20)
+       egas    = eion + eele + epos + ecoul
 
        if (.false.) then
        write(*,'(a6,e15.5)') 'egas', egas
@@ -41,7 +41,7 @@
        degasdaz = deiondaz + deepdaz + decouldaz 
        degasdzz = deiondzz + deepdzz + decouldzz 
 
-       sgas    = max(sion + sele + spos + scoul, 1d-20)
+       sgas    = sion + sele + spos + scoul
 
        dsgasdd  = dsiondd  + dsepdd  + dscouldd
        dsgasdt  = dsiondt  + dsepdt  + dscouldt

eos/private/helm.f90

@@ -183,7 +183,7 @@ subroutine ideal_eos( &
                         phase, latent_ddlnT, latent_ddlnRho, res, d_dlnd, d_dlnT, ierr)
       if(ierr/=0) return
 
-      res(i_mu) = 1d0 ! ideal assumes neutral matter, whereas pack_for_export assumes ionized matter. So we patch it up here.
+      res(i_mu) = abar ! ideal assumes neutral matter, whereas pack_for_export assumes ionized matter. So we patch it up here.
 
    end subroutine ideal_eos