Eos MESA updated tables + entropy tables for GRTDE

data/eos/mesa_opac/README

@@ -0,0 +1,10 @@
+The opac data tables are too large to be stored in the Phantom git repository
+They will be downloaded automatically when you run the code
+
+or can be retrieved manually using wget from the phantom website, e.g.:
+
+wget http://users.monash.edu.au/~dprice/phantom/data/eos/mesa_opac/opacs.bindata
+
+The files are:
+
+opacs.bindata

src/main/cons2primsolver.f90

@@ -135,7 +135,7 @@ end subroutine primitive2conservative
 subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ierr,ien_type)
  use utils_gr,     only:get_sqrtg,get_sqrt_gamma
  use metric_tools, only:unpack_metric
- use eos,          only:ieos,gmw,get_entropy,get_p_from_rho_s,gamma_global=>gamma
+ use eos,          only:ieos,gmw,get_entropy,get_p_from_rho_s,get_u_from_rho_s,gamma_global=>gamma
  use io,           only:fatal
  use physcon,      only:radconst,Rg
  use units,        only:unit_density,unit_ergg
@@ -198,6 +198,10 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie
        temp_eos = temp
        have_eos_cache = .true.
        select case(ieos)
+       case (10)
+         ! inputs and outputs are all in code units
+         call get_u_from_rho_s(ieos,en,dens,u)
+
        case (12)
           cgsdens = dens * unit_density
           cgsu = 1.5*rg*temp/gmw + radconst*temp**4/cgsdens
@@ -212,7 +216,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie
           endif
        case (2)
        case default
-          call fatal('cons2primsolver','only implemented for eos 2 and 12')
+          call fatal('cons2primsolver','only implemented for eos 2,10 and 12')
        end select
     else
        p = en*dens**gamma
@@ -261,6 +265,10 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie
        call get_p_from_rho_s(ieos,en,dens,gmw,P,temp)
     endif
     select case(ieos)
+    case (10)
+      ! inputs and outputs are all in code units
+       call get_u_from_rho_s(ieos,en,dens,u)
+
     case (12)
        cgsdens = dens * unit_density
        cgsu = 1.5*rg*temp/gmw + radconst*temp**4/cgsdens

src/main/datafiles.f90

@@ -19,6 +19,8 @@ module datafiles
 ! :Dependencies: datautils, io, mpiutils
 !
  implicit none
+! MESA EOS table version (hidden variable, set via input file)
+ integer, public :: eosmesa_version = 2   ! 0 or 1 uses the older version of the mesa tables (Reichardt et al 2020), which are no longer used by default
 
 contains
 
@@ -62,6 +64,18 @@ function map_dir_to_web(search_dir) result(url)
  !print*,' search_dir=',trim(search_dir)
  select case(search_dir)
  case('data/eos/mesa')
+    ! EOS table versions:
+    !   0,1 = legacy tables (Reichardt et al. 2020)
+    !   2   = current tables
+    select case(eosmesa_version)
+    case(0,1)
+       url = 'https://zenodo.org/records/13148447/files/'
+    case(2)
+       url = 'https://zenodo.org/records/21253459/files/'
+    case default
+       stop 'Unknown eosmesa_version'
+    end select
+ case('data/eos/mesa_opac') !!!  same link as old eos/mesa, but this is used only to download the opacity tables
     url = 'https://zenodo.org/records/13148447/files/'
  case('data/eos/shen')
     url = 'https://zenodo.org/records/13163155/files/'

src/main/eos.f90

@@ -69,7 +69,7 @@ module eos
  public  :: get_TempPresCs,get_spsound,get_temperature,get_pressure,get_cv
  public  :: eos_is_non_ideal,eos_outputs_mu,eos_outputs_gamma,eos_outputs_gasP
  public  :: eos_outputs_temp,get_local_u_internal,get_temperature_from_u
- public  :: calc_temp_and_ene,entropy,get_rho_from_p_s,get_u_from_rhoT
+ public  :: calc_temp_and_ene,entropy,get_rho_from_p_s,get_u_from_rhoT,get_u_from_rho_s
  public  :: calc_rho_from_PT,get_entropy,get_p_from_rho_s
  public  :: init_eos,finish_eos
  public  :: write_options_eos,read_options_eos,set_defaults_eos
@@ -165,9 +165,9 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam
  !
  ! Check to see if equation of state is compatible with GR cons2prim routines
  !
- if (gr .and. .not.any((/2,4,11,12/)==eos_type)) then
+ if (gr .and. .not.any((/2,4,10,11,12/)==eos_type)) then
     ponrhoi = 0.; spsoundi = 0. ! avoid compiler warning
-    call fatal('eos','GR currently only works for ieos=2,12 or 11',&
+    call fatal('eos','GR currently only works for ieos=2,4,10,12 or 11',&
          var='eos_type',val=real(eos_type))
  endif
 
@@ -1051,15 +1051,15 @@ end subroutine calc_rho_from_PT
 !-----------------------------------------------------------------------
 function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in)
  use io,                only:fatal,warning
- use physcon,           only:radconst,kb_on_mh,Rg
+ use physcon,           only:radconst,kb_on_mh,Rg, kboltz, avogadro
  use eos_idealplusrad,  only:get_idealgasplusrad_tempfrompres
  use eos_mesa,          only:get_eos_eT_from_rhop_mesa
  use mesa_microphysics, only:getvalue_mesa
  real,    intent(in) :: rho,pres,mu_in
  integer, intent(in) :: ientropy
  real,    intent(in),  optional :: eint_in,T_in,Trad_in
  integer, intent(out), optional :: ierr
- real                           :: mu,entropy,logentropy,temp,Trad,eint
+ real                           :: mu,entropy,temp,Trad,eint
 
  if (present(ierr)) ierr=0
 
@@ -1101,21 +1101,23 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in)
        call get_eos_eT_from_rhop_mesa(rho,pres,eint,temp)
     endif
 
-    ! Get entropy from rho and eint from MESA tables
+    ! Get entropy from rho and eint from MESA tables (output is not logs, it is s)
+
     if (present(ierr)) then
-       call getvalue_mesa(rho,eint,9,logentropy,ierr)
+       call getvalue_mesa(rho,eint,9,entropy,ierr)
     else
-       call getvalue_mesa(rho,eint,9,logentropy)
+       call getvalue_mesa(rho,eint,9,entropy)
     endif
-    entropy = 10.**logentropy
-
+    entropy = entropy * kboltz*avogadro ! the MESA tables are specific entropy divided by (avo*kerg).
+    ! the units of entropy with cgs inputs should be erg/g/K now
  case default
     entropy = 0.
     call fatal('eos','Unknown ientropy (can only be 1, 2, or 3)')
  end select
 
 end function entropy
 
+! input and output are in code units. entropy is in erg/g/K
 real function get_entropy(rho,pres,mu_in,ieos)
  use units,   only:unit_density,unit_pressure,unit_ergg
  use physcon, only:kboltz
@@ -1134,7 +1136,8 @@ real function get_entropy(rho,pres,mu_in,ieos)
     cgss = entropy(cgsrho,cgspres,mu_in,1)
  end select
  cgss = cgss/kboltz ! s/kb
- get_entropy = cgss/unit_ergg
+
+ get_entropy = cgss/unit_ergg ! units in erg/grK, here it turns to code units
 
 end function get_entropy
 
@@ -1170,31 +1173,38 @@ end subroutine get_rho_from_p_s
 
 !-----------------------------------------------------------------------
 !+
-!  Calculate temperature given density and entropy using Newton-Raphson
-!  method
+!  Calculate temperature and pressure given density and entropy using Newton-Raphson
+!  method (for EOS MESA it is only used in the GR case, and they are only read from the tables)
 !+
 !-----------------------------------------------------------------------
 subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out)
  use physcon, only:Rg,mass_proton_cgs
  use io,      only:fatal
- use eos_idealplusrad, only:get_idealplusrad_tempfromrhoS
+ use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_pres,&
+                            get_idealplusrad_tempfromrhoS
+ use eos_mesa,          only: get_eos_ptemp_from_rhos_mesa_gr
  use units,   only:unit_density,unit_pressure,unit_ergg
  real,    intent(in)    :: S,mu,rho
  real,    intent(inout) :: temp
  real,    intent(out)   :: P
  integer, intent(in)    :: ieos
+ real                :: cgsrho,cgspres,cgss
+ real,    parameter  :: eoserr=1e-12
+ integer, parameter  :: nitermax = 1000
  integer, intent(out), optional :: niter_out
- real :: cgsrho,cgspres,cgss
 
  ! change to cgs unit
  cgsrho = rho*unit_density
- cgss   = s*unit_ergg
+ cgss   = S*unit_ergg
  if (present(niter_out)) niter_out = 0
 
  select case (ieos)
  case (2,5)
     temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.)
     cgspres = cgsrho*Rg*temp / mu
+ case(10)
+    !!! For GR
+    call get_eos_ptemp_from_rhos_mesa_gr(cgsrho,cgss,cgspres,temp)
  case (12)
     call get_idealplusrad_tempfromrhoS(cgsrho,cgss,mu,temp,cgspres,niter_out)
  case default
@@ -1211,6 +1221,39 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out)
 
 end subroutine get_p_from_rho_s
 
+!-----------------------------------------------------------------------
+!+
+!  Calculate temperature given density and entropy using EOS MESA tables for GR case
+!+
+!-----------------------------------------------------------------------
+subroutine get_u_from_rho_s(ieos,S,rho,u)
+ use io,      only:fatal
+ use units,   only:unit_density,unit_ergg
+ use eos_mesa,          only: get_eos_u_from_rhos_mesa_gr
+ real,    intent(in)    :: S,rho
+ real,    intent(out)   :: u
+ integer, intent(in)    :: ieos
+ real                :: cgsrho,cgss, cgsu
+
+ ! change to cgs unit
+ cgsrho = rho*unit_density
+ cgss   = S*unit_ergg
+
+ select case (ieos)
+ case(10)
+    !!! For GR
+    call get_eos_u_from_rhos_mesa_gr(cgsrho,cgss,cgsu)
+
+ case default
+    cgsu = 0.
+    call fatal('eos','[get_u_from_rho_s] only implemented for eos 10')
+ end select
+
+ ! change back to code unit
+ u = cgsu / unit_ergg
+
+end subroutine get_u_from_rho_s
+
 !-----------------------------------------------------------------------
 !+
 !  Calculate mean molecular weight from X and Z, assuming complete
@@ -1859,7 +1902,6 @@ subroutine set_defaults_eos
  C_ent              = 3.
  polyk2             = 0. ! only used for ieos=8
  use_var_comp = .false.  ! variable composition
-
 end subroutine set_defaults_eos
 
 end module eos

src/main/eos_mesa.f90

@@ -14,11 +14,11 @@ module eos_mesa
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: mesa_microphysics, physcon
+! :Dependencies: mesa_microphysics, physcon, dim
 !
 
  use mesa_microphysics
-
+ use dim,           only:gr
  implicit none
  logical, private :: mesa_initialised = .false.
 
@@ -58,7 +58,18 @@ subroutine init_eos_mesa(x,z,ierr)
  call get_eos_constants_mesa(ierr)
  if (ierr /= 0) return
 
+ !!! only read GR tables if it is a GR run
+ mesa_eos_gr_prefix="output_rhos_"
+ if (gr) then
+    call get_eos_constants_mesa_gr(ierr)
+    if (ierr /= 0) return
+
+    call read_eos_mesa_gr(x,z,ierr)
+    if (ierr /= 0) return
+ end if
+
  call read_eos_mesa(x,z,ierr)
+ if (ierr /= 0) return
  call get_opacity_constants_mesa
  call read_opacity_mesa(x,z)
 
@@ -134,6 +145,37 @@ subroutine get_eos_pressure_temp_mesa(den,eint,pres,temp)
 
 end subroutine get_eos_pressure_temp_mesa
 
+!----------------------------------------------------------------
+!+
+!  subroutine returns pressure and temperature as
+!  a function of density/entropy for GR tables
+!+
+!----------------------------------------------------------------
+subroutine get_eos_ptemp_from_rhos_mesa_gr(den,s,pres,temp)
+ real, intent(in)  :: den, s
+ real, intent(out) :: pres, temp
+
+ call getvalue_mesa_gr(den,s,1,pres)
+ call getvalue_mesa_gr(den,s,3,temp)
+
+end subroutine get_eos_ptemp_from_rhos_mesa_gr
+
+!----------------------------------------------------------------
+!+
+!  subroutine returns internal energy as
+!  a function of density/entropy for GR tables
+!+
+!----------------------------------------------------------------
+subroutine get_eos_u_from_rhos_mesa_gr(den,s,u)
+ real, intent(in)  :: den, s
+ real, intent(out) :: u
+
+ call getvalue_mesa_gr(den,s,2,u)
+
+end subroutine get_eos_u_from_rhos_mesa_gr
+
+
+
 !----------------------------------------------------------------
 !+
 !  subroutine returns internal energy and temperature from