diff --git a/data/eos/mesa_opac/README b/data/eos/mesa_opac/README new file mode 100644 index 000000000..241cd9001 --- /dev/null +++ b/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 diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 7a616d1f7..65d7c1569 100644 --- a/src/main/cons2primsolver.f90 +++ b/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 diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index 7cfdd0169..efad21ba2 100644 --- a/src/main/datafiles.f90 +++ b/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/' diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 01c4da0a0..75d090d6c 100644 --- a/src/main/eos.f90 +++ b/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,7 +1051,7 @@ 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 @@ -1059,7 +1059,7 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_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,14 +1101,15 @@ 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)') @@ -1116,6 +1117,7 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in) 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 diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index 985996846..066ac80eb 100644 --- a/src/main/eos_mesa.f90 +++ b/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 diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 827b9e647..a2931c9bc 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -42,6 +42,18 @@ module mesa_microphysics real, allocatable :: mesa_eos0(:,:,:,:) real, allocatable :: mesa_de_data0(:,:,:,:,:) +! EOS tables in entropy and density for GR + integer :: mesa_eos_gr_nz=3, mesa_eos_gr_nh=5 !Number of Z values in EOS tables, number of hydrogen values in EOS tables + real :: mesa_eos_gr_z1=0.d0, mesa_eos_gr_h1=0.d0, mesa_eos_gr_dz=0.02d0, mesa_eos_gr_dh=0.2d0 !lowest Z, lowest H, delta Z, delta H in EOS tables + integer :: mesa_eos_gr_ns, mesa_eos_gr_nrho, mesa_eos_gr_nvar2 + real :: mesa_eos_gr_rho1, mesa_eos_gr_s1, mesa_eos_gr_ds, mesa_eos_gr_drho + character(len=200) :: mesa_eos_gr_prefix="eos_gr_" ! prefix for EOS GR table files + real, allocatable :: mesa_eos_gr_z(:), mesa_eos_gr_h(:),mesa_eos_gr_logSs(:), mesa_eos_gr_logrhos(:) + integer, allocatable :: mesa_eos_gr_data_exists(:,:) + real, allocatable :: mesa_gr_ds_data(:,:,:) + real, allocatable :: mesa_eos_gr0(:,:,:,:) + real, allocatable :: mesa_gr_ds_data0(:,:,:,:,:) + public :: get_opacity_constants_mesa public :: read_opacity_mesa public :: get_kappa_mesa @@ -53,6 +65,8 @@ module mesa_microphysics !public :: getgamma1_mesa !public :: getpressure_mesa public :: getvalue_mesa + public :: read_eos_mesa_gr + public :: getvalue_mesa_gr private :: eos_cubic_spline_mesa private :: cubic_spline_mesa @@ -66,7 +80,7 @@ subroutine get_opacity_constants_mesa ! Find the opacity table filename = 'opacs'//trim(mesa_opacs_suffix)//'.bindata' - opacs_file = find_phantom_datafile(filename,'eos/mesa') + opacs_file = find_phantom_datafile(filename,'eos/mesa_opac') ! Read the constants from the header of the opacity file open(newunit=fnum,file=trim(opacs_file),status='old',action='read',form='unformatted') @@ -100,7 +114,7 @@ subroutine read_opacity_mesa(x,z) !Find the opacity table filename = trim(mesa_opacs_dir)//'opacs'//trim(mesa_opacs_suffix)//'.bindata' ! filename = trim(mesa_opacs_dir)//'/'//'opacs'//trim(mesa_opacs_suffix)//'.bindata' - opacs_file = find_phantom_datafile(filename,'eos/mesa') + opacs_file = find_phantom_datafile(filename,'eos/mesa_opac') open(unit=fnum,file=trim(opacs_file),status='old',action='read',form='unformatted') read(fnum) mesa_opacs_nz,mesa_opacs_nx,mesa_opacs_nr,mesa_opacs_nt @@ -312,8 +326,64 @@ subroutine get_eos_constants_mesa(ierr) allocate(mesa_eos0(mesa_eos_nz,mesa_eos_nh,mesa_eos_ne,mesa_eos_nv)) allocate(mesa_de_data0(mesa_eos_nz,mesa_eos_nh,mesa_eos_ne,mesa_eos_nv,mesa_eos_nvar2)) + return + end subroutine get_eos_constants_mesa + +! Get the constants to be used in the MESA EoS for GR +subroutine get_eos_constants_mesa_gr(ierr) + integer, intent(out) :: ierr + character(len=20) :: zz, hh + character(len=300) :: filename + character(len=8) :: fmt1 + integer :: fnum, i + + ! Allocate the arrays to carry X and Z values. + allocate(mesa_eos_gr_z(mesa_eos_gr_nz), mesa_eos_gr_h(mesa_eos_gr_nh)) + + ! Assigning the X and Z values of the EoS tables to be read in (based on parameters at the beginning of the file. + mesa_eos_gr_z(1)=mesa_eos_gr_z1 + mesa_eos_gr_h(1)=mesa_eos_gr_h1 + + do i=2,mesa_eos_gr_nz + mesa_eos_gr_z(i)=mesa_eos_gr_z(i-1)+mesa_eos_gr_dz + enddo + + do i=2,mesa_eos_gr_nh + mesa_eos_gr_h(i)=mesa_eos_gr_h(i-1)+mesa_eos_gr_dh + enddo + + fmt1 = '(F4.2)' ! + + fnum=126 + + ! Convert first X and Z to string + write (zz,fmt1) mesa_eos_gr_z1 + write (hh,fmt1) mesa_eos_gr_h1 + + ! Find the first EoS tables + filename = trim(mesa_eos_gr_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' + filename = find_phantom_datafile(filename,'eos/mesa') + + ! Read constants from the header of first EoS tables + open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted',iostat=ierr) +! allocate GR tables + if (ierr /= 0) return + read(fnum) mesa_eos_gr_ns, mesa_eos_gr_nrho, mesa_eos_gr_nvar2 + close(fnum) + + allocate(mesa_eos_gr_logSs(mesa_eos_gr_ns), mesa_eos_gr_logrhos(mesa_eos_gr_nrho)) + allocate(mesa_eos_gr_data_exists(mesa_eos_gr_nz,mesa_eos_gr_nh)) + allocate(mesa_gr_ds_data(mesa_eos_gr_ns,mesa_eos_gr_nrho,mesa_eos_gr_nvar2)) + allocate(mesa_eos_gr0(mesa_eos_gr_nz,mesa_eos_gr_nh,mesa_eos_gr_ns,mesa_eos_gr_nrho)) + allocate(mesa_gr_ds_data0(mesa_eos_gr_nz,mesa_eos_gr_nh,mesa_eos_gr_ns,mesa_eos_gr_nrho,mesa_eos_gr_nvar2)) + + return + +end subroutine get_eos_constants_mesa_gr + + ! Read MESA EoS tables, and then construct a new array for the specific values of X and Z subroutine read_eos_mesa(x,z,ierr) real, intent(in) :: x, z @@ -347,13 +417,13 @@ subroutine read_eos_mesa(x,z,ierr) ! Find the EoS tables filename = trim(mesa_eos_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' filename = find_phantom_datafile(filename,'eos/mesa') - ! Note that the data exists mesa_eos_data_exists(i,j)=1 ! Read in the size of the tables and the data ! i and j hold the Z and X values respectively ! k, l and m hold the values of V, Eint and the data respectively + open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted') read(fnum) mesa_eos_ne, mesa_eos_nv, mesa_eos_nvar2 read(fnum)(mesa_eos_logVs(k),k=1,mesa_eos_nv) @@ -419,12 +489,117 @@ subroutine read_eos_mesa(x,z,ierr) end subroutine read_eos_mesa +! same as above but for GR tables (entropy and density instead of internal energy and density) +subroutine read_eos_mesa_gr(x,z,ierr) + real, intent(in) :: x, z + integer, intent(out) :: ierr + real, parameter :: arad=7.5657d-15 + integer :: i,j,k,l,m + character(len=300) :: filename + character(len=8) :: fmt1 + character(len=20) :: zz, hh + integer :: nz1,nz2,nx1,nx2 + real :: dz,dx + integer :: fnum + + fmt1 = '(F4.2)' + + fnum = 126 + ierr = 0 + + ! following lines to prevent compiler warnings + nx1 = 0 + nx2 = 0 + dx = 0. + + !Loop over files + do i=1,mesa_eos_gr_nz + write (zz,fmt1) mesa_eos_gr_z(i) ! Z value as string + + do j=1,mesa_eos_gr_nh + write (hh,fmt1) mesa_eos_gr_h(j) ! X value as string + ! Find the EoS tables + filename = trim(mesa_eos_gr_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' + filename = find_phantom_datafile(filename,'eos/mesa') + + ! Note that the data exists + mesa_eos_gr_data_exists(i,j)=1 + + ! Read in the size of the tables and the data + ! i and j hold the Z and X values respectively + ! k, l and m hold the values of V, Eint and the data respectively + open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted') + read(fnum) mesa_eos_gr_ns, mesa_eos_gr_nrho, mesa_eos_gr_nvar2 + read(fnum)(mesa_eos_gr_logrhos(k),k=1,mesa_eos_gr_nrho) + read(fnum)(mesa_eos_gr_logss(l),l=1,mesa_eos_gr_ns) + do k=1,mesa_eos_gr_nrho + do l=1,mesa_eos_gr_ns + read(fnum) (mesa_gr_ds_data0(i,j,l,k,m),m=1,mesa_eos_gr_nvar2) + enddo + enddo + + close(fnum) + enddo + enddo + + ! Chooses the Z values from the available values which is greater than the input Z + nz2=mesa_eos_gr_nz + do i=2,mesa_eos_gr_nz-1 + if (mesa_eos_gr_z(i) >= z) then + nz2=i + exit + endif + enddo + + ! Chooses Z values below input Z, and then finds fraction of input Z distance between two enclosing Z (sorry if confusing, just read the code) + nz1 = nz2 - 1 + dz=(z-mesa_eos_gr_z(nz1))/(mesa_eos_gr_z(nz2)-mesa_eos_gr_z(nz1)) + + ! Seems to do the same as above, except for X. I can't quite figure out the weird if statements. + do j=2,mesa_eos_gr_nh + if (mesa_eos_gr_h(j) > x.or.j==mesa_eos_gr_nh) then + nx2=j + if (j==2.and.mesa_eos_gr_data_exists(i,j-2)==0) then + nx2=3 + elseif (j==mesa_eos_gr_nh.and.mesa_eos_gr_data_exists(i,mesa_eos_gr_nh-1)==0) then + nx2=mesa_eos_gr_nh-1 + endif + exit + endif + enddo + nx1 = nx2 - 1 + dx = (x-mesa_eos_gr_h(nx1))/(mesa_eos_gr_h(nx2)-mesa_eos_gr_h(nx1)) + + ! Bilinear interpolation of values from the four tables (two values each of X and Z) + do l=1,mesa_eos_gr_ns + do k=1,mesa_eos_gr_nrho + do m=1,mesa_eos_gr_nvar2 + mesa_gr_ds_data(l,k,m)=(1.d0-dx)*(1.d0-dz)*mesa_gr_ds_data0(nz1,nx1,l,k,m)+dz*(1.d0-dx)*mesa_gr_ds_data0(nz2,nx1,l,k,m) & + +dx*(1.d0-dz)*mesa_gr_ds_data0(nz1,nx2,l,k,m)+dx*dz*mesa_gr_ds_data0(nz2,nx2,l,k,m) + enddo + enddo + enddo + + ! Save some things to make the interpolation fast + mesa_eos_gr_drho=mesa_eos_gr_logrhos(2)-mesa_eos_gr_logrhos(1) + mesa_eos_gr_ds=mesa_eos_gr_logss(2)-mesa_eos_gr_logss(1) + mesa_eos_gr_rho1=mesa_eos_gr_logrhos(1) + mesa_eos_gr_s1=mesa_eos_gr_logss(1) + + !deallocate(mesa_eos_gr0) + !deallocate(mesa_gr_ds_data0) + + return + +end subroutine read_eos_mesa_gr + ! Get value from the MESA tables given a specific value of density and internal energy ! The columns in the data are: ! 1. logRho 2. logP 3. logPgas 4. logT ! 5. dlnP/dlnrho|e 6. dlnP/dlne|rho 7. dlnT/dlnrho|e 8. dlnT/dlne|rho ! 9. logS 10. dlnT/dlnP|S 11. Gamma1 12. gamma ! Note: ivout=1,2,3,4 returns the unlogged quantity + pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) real, intent(in) :: rho, eint real, intent(out) :: vout @@ -454,16 +629,18 @@ pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) dv = (logv - mesa_eos_logVs(nv)) / mesa_eos_dv ! If the given Eint and V fall within the limits of the table, then use cubic spline interpolation to find value - ! Else use linear extrapolation beyond the limits of the table (I think? I can't tell if this is the correct way to do this) + ! Else use linear extrapolation beyond the limits of the table if (ne > 1 .and. nv > 1 .and. ne < mesa_eos_ne-1 .and. nv < mesa_eos_nv-1) then call eos_cubic_spline_mesa(ne,nv,loge,logv,ivout,vout,nx,dx) - if (ivout < 5) vout=10.d0**vout + if (ivout < 5 .or. ivout == 9) vout=10.d0**vout ! this is applied only for the first four columns of the data which are stored in log10 form. if (present(ierr)) ierr = 0 else - vout = 10.d0**((1.d0-de) * (1.d0-dv) * mesa_de_data(ne,nv,ivout) + & - de * (1.d0-dv) * mesa_de_data(ne+1,nv,ivout) + & - (1.d0-de) * dv * mesa_de_data(ne,nv+1,ivout) + & - de * dv * mesa_de_data(ne+1,nv+1,ivout)) + vout = (1.d0-de) * (1.d0-dv) * mesa_de_data(ne,nv,ivout) + & + de * (1.d0-dv) * mesa_de_data(ne+1,nv,ivout) + & + (1.d0-de) * dv * mesa_de_data(ne,nv+1,ivout) + & + de * dv * mesa_de_data(ne+1,nv+1,ivout) + + if (ivout <= 4 .or. ivout == 9) vout = 10.d0**vout if (present(ierr)) ierr = 1 ! warn if extrapolating endif @@ -471,6 +648,58 @@ pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) end subroutine getvalue_mesa +! Get value from the MESA tables for GR given a specific value of density and entropy +! The columns in the data are: +! 1. logP 2. logu 3. logT 4. Gamma1 +! Note: ivout=1,2,3,4 returns the unlogged quantity +! inputted s should be in cgs units, then it has to be converted to log10(S/(k_B*N_A)) in cgs units (because the MESA tables are in these units) +pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) + use physcon, only: kboltz,avogadro + real, intent(in) :: rho, s + real, intent(out) :: vout + integer, intent(in) :: ivout + integer, intent(out), optional :: ierr + real :: logs, logrho, ds, drho + integer :: ns, nrho + real :: dx + integer :: nx + + logs = log10(s/(kboltz*avogadro)) ! convert to log10(S/(k_B*N_A)) in cgs units (because the MESA tables are in these units) + logrho = log10(rho) + + ! Get the s and rho indices for looking up the tables + ns = 1 + int((logs - mesa_eos_gr_s1) / mesa_eos_gr_ds) + nrho = 1 + int((logrho - mesa_eos_gr_rho1) / mesa_eos_gr_drho) + + ! Allow extrapolation + if (ns < 1) ns=1 + if (ns > mesa_eos_gr_ns-1) ns=mesa_eos_gr_ns-1 + if (nrho < 1) nrho=1 + if (nrho > mesa_eos_gr_nrho-1) nrho=mesa_eos_gr_nrho-1 + + ds = (logs - mesa_eos_gr_logSs(ns)) / mesa_eos_gr_ds + drho = (logrho - mesa_eos_gr_logRhos(nrho)) / mesa_eos_gr_drho + + ! If the given S and Rho fall within the limits of the table, then use cubic spline interpolation to find value + ! Else use linear extrapolation beyond the limits of the table + if (ns > 1 .and. nrho > 1 .and. ns < mesa_eos_gr_ns-1 .and. nrho < mesa_eos_gr_nrho-1) then + call eos_cubic_spline_mesa_gr(ns,nrho,logs,logrho,ivout,vout,nx,dx) + if (ivout < 4) vout=10.d0**vout ! this is applied only for the first three columns of the data which are stored in log10 form. + if (present(ierr)) ierr = 0 + else + vout = (1.d0-ds) * (1.d0-drho) * mesa_gr_ds_data(ns,nrho,ivout) + & + ds * (1.d0-drho) * mesa_gr_ds_data(ns+1,nrho,ivout) + & + (1.d0-ds) * drho * mesa_gr_ds_data(ns,nrho+1,ivout) + & + ds * drho * mesa_gr_ds_data(ns+1,nrho+1,ivout) + if (ivout <= 3) vout = 10.d0**vout + if (present(ierr)) ierr = 1 ! warn if extrapolating + endif + + return + +end subroutine getvalue_mesa_gr + + !only use if between e(2) < e < e(n_e-1) and v(2) < v < v(n_v-1) pure subroutine eos_cubic_spline_mesa(e1,v1,e,v,n_var,z,h1,dh) @@ -527,6 +756,64 @@ pure subroutine eos_cubic_spline_mesa(e1,v1,e,v,n_var,z,h1,dh) end subroutine eos_cubic_spline_mesa +!Spline for GR: since eos_cubic_spline_mesa has the original tables embedded in it, I have to make a new one for the GR tables. +! It's basically the same as above but with different variables and using mesa_gr_ds_data instead of mesa_de_data +! only use if between s(2) < s < s(n_s-1) and rho(2) < rho < rho(n_rho-1) + +pure subroutine eos_cubic_spline_mesa_gr(s1,rho1,s,rho,n_var,z,h1,dh) + +! use mesa_eos_gr_logss, mesa_eos_gr_logrhos + implicit none + + integer, intent(in) :: s1, rho1, h1, n_var + real, intent(in) :: s, rho, dh + + real, intent(out) :: z + + real :: x0,x1,x2,x3,y0,y1,y2,y3 + real :: as,bs,tt + + real :: yy(4) + integer :: j + + x0=mesa_eos_gr_logss(s1-1) + x1=mesa_eos_gr_logss(s1+0) + x2=mesa_eos_gr_logss(s1+1) + x3=mesa_eos_gr_logss(s1+2) + + tt=(s-x1)/mesa_eos_gr_ds + + do j=1,4 + + y0 = mesa_gr_ds_data(s1-1,rho1+j-2,n_var) + y1 = mesa_gr_ds_data(s1+0,rho1+j-2,n_var) + y2 = mesa_gr_ds_data(s1+1,rho1+j-2,n_var) + y3 = mesa_gr_ds_data(s1+2,rho1+j-2,n_var) + + call cubic_spline_mesa(x0,x1,x2,x3,y0,y1,y2,y3,as,bs) + yy(j)=(1.d0-tt)*y1+tt*y2+tt*(1.d0-tt)*(as*(1.d0-tt)+bs*tt) + + enddo + + x0=mesa_eos_gr_logrhos(rho1-1) + x1=mesa_eos_gr_logrhos(rho1+0) + x2=mesa_eos_gr_logrhos(rho1+1) + x3=mesa_eos_gr_logrhos(rho1+2) + + y0=yy(1) + y1=yy(2) + y2=yy(3) + y3=yy(4) + + tt=(rho-x1)/mesa_eos_gr_drho + + call cubic_spline_mesa(x0,x1,x2,x3,y0,y1,y2,y3,as,bs) + z=(1.d0-tt)*y1+tt*y2+tt*(1.d0-tt)*(as*(1.d0-tt)+bs*tt) + + return + +end subroutine eos_cubic_spline_mesa_gr + pure subroutine cubic_spline_mesa(x0,x1,x2,x3,y0,y1,y2,y3,as,bs) implicit none @@ -563,7 +850,18 @@ subroutine deallocate_arrays_mesa if (allocated(mesa_eos0)) deallocate(mesa_eos0) if (allocated(mesa_de_data)) deallocate(mesa_de_data) if (allocated(mesa_de_data0)) deallocate(mesa_de_data0) +!eos gr + if (allocated(mesa_eos_gr_z)) deallocate(mesa_eos_gr_z) + if (allocated(mesa_eos_gr_h)) deallocate(mesa_eos_gr_h) + if (allocated(mesa_eos_gr0)) deallocate(mesa_eos_gr0) + if (allocated(mesa_eos_gr_logss)) deallocate(mesa_eos_gr_logss) + if (allocated(mesa_eos_gr_logrhos)) deallocate(mesa_eos_gr_logrhos) + if (allocated(mesa_eos_gr_data_exists)) deallocate(mesa_eos_gr_data_exists) + if (allocated(mesa_gr_ds_data)) deallocate(mesa_gr_ds_data) + if (allocated(mesa_gr_ds_data0)) deallocate(mesa_gr_ds_data0) + end subroutine deallocate_arrays_mesa + end module mesa_microphysics