From bd792a65c09ab23f9d67640e6af23f97a65d77e1 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Thu, 11 Jun 2026 12:16:03 +1000 Subject: [PATCH 01/42] I started adding modules to import eos mesa tables in entropy and density for GR, I want to change as little as possible though and preferably get away with GR users setting the same EOS number too --- src/main/eos_mesa_microphysics.f90 | 116 +++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 827b9e647..c099e4aec 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -42,6 +42,15 @@ module mesa_microphysics real, allocatable :: mesa_eos0(:,:,:,:) real, allocatable :: mesa_de_data0(:,:,:,:,:) +! EOS tables in entropy and density for GR + integer :: mesa_eos_ns, mesa_eos_nrho !, mesa_eos_nvar2 + real :: mesa_eos_rho1, mesa_eos_s1, mesa_eos_ds, mesa_eos_drho + real, allocatable :: mesa_eos_logss(:), mesa_eos_logrhos(:) +! integer, allocatable :: mesa_eos_data_exists(:,:) + real, allocatable :: mesa_ds_data(:,:,:) +! real, allocatable :: mesa_eos0(:,:,:,:) + real, allocatable :: mesa_ds_data0(:,:,:,:,:) + public :: get_opacity_constants_mesa public :: read_opacity_mesa public :: get_kappa_mesa @@ -419,12 +428,119 @@ 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 = 120 + ierr = 0 + + ! following lines to prevent compiler warnings + nx1 = 0 + nx2 = 0 + dx = 0. + + !Loop over files + do i=1,mesa_eos_nz + write (zz,fmt1) mesa_eos_z(i) ! Z value as string + + do j=1,mesa_eos_nh + write (hh,fmt1) mesa_eos_h(j) ! X value as string +!!!!!!!!!! have to rename these variables to not conflict with the above subroutine, but I can't be bothered to do that rn, so just ignore the fact that they have the same names as the above subroutine !!!!!!!!!! + ! 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_ns, mesa_eos_nrho, mesa_eos_nvar2 + read(fnum)(mesa_eos_logrhos(k),k=1,mesa_eos_nrho) + read(fnum)(mesa_eos_logss(l),l=1,mesa_eos_ns) + do k=1,mesa_eos_nrho + do l=1,mesa_eos_ns + read(fnum) (mesa_ds_data0(i,j,l,k,m),m=1,mesa_eos_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_nz + do i=2,mesa_eos_nz-1 + if (mesa_eos_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_z(nz1))/(mesa_eos_z(nz2)-mesa_eos_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_nh + if (mesa_eos_h(j) > x.or.j==mesa_eos_nh) then + nx2=j + if (j==2.and.mesa_eos_data_exists(i,j-2)==0) then + nx2=3 + elseif (j==mesa_eos_nh.and.mesa_eos_data_exists(i,mesa_eos_nh-1)==0) then + nx2=mesa_eos_nh-1 + endif + exit + endif + enddo + nx1 = nx2 - 1 + dx = (x-mesa_eos_h(nx1))/(mesa_eos_h(nx2)-mesa_eos_h(nx1)) + + ! Bilinear interpolation of values from the four tables (two values each of X and Z) + do l=1,mesa_eos_ns + do k=1,mesa_eos_nrho + do m=1,mesa_eos_nvar2 + mesa_ds_data(l,k,m)=(1.d0-dx)*(1.d0-dz)*mesa_ds_data0(nz1,nx1,l,k,m)+dz*(1.d0-dx)*mesa_ds_data0(nz2,nx1,l,k,m) & + +dx*(1.d0-dz)*mesa_ds_data0(nz1,nx2,l,k,m)+dx*dz*mesa_ds_data0(nz2,nx2,l,k,m) + enddo + enddo + enddo + + ! Save some things to make the interpolation fast + mesa_eos_drho=mesa_eos_logrhos(2)-mesa_eos_logrhos(1) + mesa_eos_ds=mesa_eos_logss(2)-mesa_eos_logss(1) + mesa_eos_rho1=mesa_eos_logrhos(1) + mesa_eos_s1=mesa_eos_logss(1) + + !deallocate(mesa_eos0) + !deallocate(mesa_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 + +!!! I don't think a GR version of this is needed, the interpolation is the same. pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) real, intent(in) :: rho, eint real, intent(out) :: vout From b2be5ff4a2798e72231da21df52a4c21eda01e7c Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Thu, 11 Jun 2026 14:23:11 +1000 Subject: [PATCH 02/42] adding some skeletons of some functions needed, just saving now --- src/main/eos.f90 | 42 ++++++++++++++++++++++++++++++ src/main/eos_mesa.f90 | 2 ++ src/main/eos_mesa_microphysics.f90 | 38 ++++++++++++++++++++------- 3 files changed, 73 insertions(+), 9 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index d75ad799b..8801f8374 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1195,6 +1195,8 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) case (2,5) temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.) cgspres = cgsrho*Rg*temp / mu + case(10) + !!! Ali add things for GR!!! case (12) call get_idealplusrad_tempfromrhoS(cgsrho,cgss,mu,temp,cgspres,niter_out) case default @@ -1211,6 +1213,46 @@ 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 Newton-Raphson +! method +!+ +!----------------------------------------------------------------------- +subroutine get_u_from_rho_s(ieos,S,rho,mu,u) + use physcon, only:radconst,Rg,mass_proton_cgs,kboltz + use io, only:fatal + use units, only:unit_density,unit_pressure,unit_ergg + real, intent(in) :: S,mu,rho + real, intent(out) :: u + integer, intent(in) :: ieos + real :: corr,df,f,cgsrho,cgsp,cgss + real, parameter :: eoserr=1e-12 + integer :: niter + integer, parameter :: nitermax = 1000 + + ! change to cgs unit + cgsrho = rho*unit_density + cgss = s*unit_ergg + + niter = 0 + select case (ieos) + case(10) + !!! Ali add things for GR!!! + + case default + cgsP = 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 diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index 985996846..3632de54f 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -61,6 +61,8 @@ subroutine init_eos_mesa(x,z,ierr) call read_eos_mesa(x,z,ierr) call get_opacity_constants_mesa call read_opacity_mesa(x,z) + mesa_eos_prefix="output_rhos_" + call read_eos_gr_mesa(x,z,ierr) ! read GR tables end subroutine init_eos_mesa diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index c099e4aec..c3ba9fbbd 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -43,10 +43,11 @@ module mesa_microphysics real, allocatable :: mesa_de_data0(:,:,:,:,:) ! EOS tables in entropy and density for GR - integer :: mesa_eos_ns, mesa_eos_nrho !, mesa_eos_nvar2 + integer :: mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 real :: mesa_eos_rho1, mesa_eos_s1, mesa_eos_ds, mesa_eos_drho + character(len=200) :: mesa_eos_gr_prefix="eos_gr_" ! prefix for EOS GR table files real, allocatable :: mesa_eos_logss(:), mesa_eos_logrhos(:) -! integer, allocatable :: mesa_eos_data_exists(:,:) + integer, allocatable :: mesa_eos_gr_data_exists(:,:) real, allocatable :: mesa_ds_data(:,:,:) ! real, allocatable :: mesa_eos0(:,:,:,:) real, allocatable :: mesa_ds_data0(:,:,:,:,:) @@ -320,6 +321,17 @@ subroutine get_eos_constants_mesa(ierr) allocate(mesa_de_data(mesa_eos_ne,mesa_eos_nv,mesa_eos_nvar2)) 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)) +! allocate GR tables (again, I wonder if it's best to make another subroutine for this) + if (ierr /= 0) return + read(fnum) mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 + close(fnum) + + allocate(mesa_eos_logss(mesa_eos_ns), mesa_eos_logrhos(mesa_eos_nrho)) + allocate(mesa_eos_gr_data_exists(mesa_eos_nz,mesa_eos_nh)) + allocate(mesa_ds_data(mesa_eos_ns,mesa_eos_nrho,mesa_eos_gr_nvar2)) + allocate(mesa_ds_data0(mesa_eos_nz,mesa_eos_nh,mesa_eos_ns,mesa_eos_nrho,mesa_eos_gr_nvar2)) + + return end subroutine get_eos_constants_mesa @@ -459,22 +471,22 @@ subroutine read_eos_mesa_gr(x,z,ierr) write (hh,fmt1) mesa_eos_h(j) ! X value as string !!!!!!!!!! have to rename these variables to not conflict with the above subroutine, but I can't be bothered to do that rn, so just ignore the fact that they have the same names as the above subroutine !!!!!!!!!! ! Find the EoS tables - filename = trim(mesa_eos_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' + 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_data_exists(i,j)=1 + 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_ns, mesa_eos_nrho, mesa_eos_nvar2 + read(fnum) mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 read(fnum)(mesa_eos_logrhos(k),k=1,mesa_eos_nrho) read(fnum)(mesa_eos_logss(l),l=1,mesa_eos_ns) do k=1,mesa_eos_nrho do l=1,mesa_eos_ns - read(fnum) (mesa_ds_data0(i,j,l,k,m),m=1,mesa_eos_nvar2) + read(fnum) (mesa_ds_data0(i,j,l,k,m),m=1,mesa_eos_gr_nvar2) enddo enddo @@ -499,9 +511,9 @@ subroutine read_eos_mesa_gr(x,z,ierr) do j=2,mesa_eos_nh if (mesa_eos_h(j) > x.or.j==mesa_eos_nh) then nx2=j - if (j==2.and.mesa_eos_data_exists(i,j-2)==0) then + if (j==2.and.mesa_eos_gr_data_exists(i,j-2)==0) then nx2=3 - elseif (j==mesa_eos_nh.and.mesa_eos_data_exists(i,mesa_eos_nh-1)==0) then + elseif (j==mesa_eos_nh.and.mesa_eos_gr_data_exists(i,mesa_eos_nh-1)==0) then nx2=mesa_eos_nh-1 endif exit @@ -513,7 +525,7 @@ subroutine read_eos_mesa_gr(x,z,ierr) ! Bilinear interpolation of values from the four tables (two values each of X and Z) do l=1,mesa_eos_ns do k=1,mesa_eos_nrho - do m=1,mesa_eos_nvar2 + do m=1,mesa_eos_gr_nvar2 mesa_ds_data(l,k,m)=(1.d0-dx)*(1.d0-dz)*mesa_ds_data0(nz1,nx1,l,k,m)+dz*(1.d0-dx)*mesa_ds_data0(nz2,nx1,l,k,m) & +dx*(1.d0-dz)*mesa_ds_data0(nz1,nx2,l,k,m)+dx*dz*mesa_ds_data0(nz2,nx2,l,k,m) enddo @@ -679,7 +691,15 @@ 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 (I wonder if it's best to make another subroutine for deallocating the GR tables) + if (allocated(mesa_eos_logss)) deallocate(mesa_eos_logss) + if (allocated(mesa_eos_logrhos)) deallocate(mesa_eos_logrhos) + if (allocated(mesa_eos_gr_data_exists)) deallocate(mesa_eos_gr_data_exists) + if (allocated(mesa_ds_data)) deallocate(mesa_ds_data) + if (allocated(mesa_ds_data0)) deallocate(mesa_ds_data0) + end subroutine deallocate_arrays_mesa + end module mesa_microphysics From 1d10bd999cc3454939c1cb9fa88cd0a062a244d3 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Thu, 11 Jun 2026 17:11:30 +1000 Subject: [PATCH 03/42] all the subroutines necessary for entropy tables to be read are added now, in eos, eos_mesa and eos_mesa_microphysics, now I need to add the entropy tables here, and also see if they get allocated correctly, I am only adding 4 variables in these tables so they occupy less size, other than that the compositions will be the same, but the minimum density should go as low as 1e-25 g/cm as fitz' simulations show. also this is not tested so I should do that --- src/main/eos.f90 | 27 ++++--- src/main/eos_mesa.f90 | 31 ++++++++ src/main/eos_mesa_microphysics.f90 | 111 ++++++++++++++++++++++++++++- 3 files changed, 153 insertions(+), 16 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 8801f8374..ce08d84d3 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1170,8 +1170,8 @@ 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) @@ -1196,7 +1196,9 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.) cgspres = cgsrho*Rg*temp / mu case(10) - !!! Ali add things for GR!!! + !!! For GR, not tested yet! + call get_eos_ptemp_from_rhos_mesa_gr(cgsrho,cgss,cgsP,temp) + case (12) call get_idealplusrad_tempfromrhoS(cgsrho,cgss,mu,temp,cgspres,niter_out) case default @@ -1215,33 +1217,28 @@ end subroutine get_p_from_rho_s !----------------------------------------------------------------------- !+ -! Calculate temperature given density and entropy using Newton-Raphson -! method +! Calculate temperature given density and entropy using EOS MESA tables for GR case !+ !----------------------------------------------------------------------- -subroutine get_u_from_rho_s(ieos,S,rho,mu,u) - use physcon, only:radconst,Rg,mass_proton_cgs,kboltz +subroutine get_u_from_rho_s(ieos,S,rho,u) use io, only:fatal use units, only:unit_density,unit_pressure,unit_ergg - real, intent(in) :: S,mu,rho + real, intent(in) :: S,rho real, intent(out) :: u integer, intent(in) :: ieos - real :: corr,df,f,cgsrho,cgsp,cgss - real, parameter :: eoserr=1e-12 - integer :: niter - integer, parameter :: nitermax = 1000 + real :: cgsrho,cgsp,cgss ! change to cgs unit cgsrho = rho*unit_density cgss = s*unit_ergg - niter = 0 select case (ieos) case(10) - !!! Ali add things for GR!!! + !!! For GR, not tested yet! + call get_eos_u_from_rhos_mesa_gr(cgsrho,cgss,cgsu) case default - cgsP = 0. + cgsu = 0. call fatal('eos','[get_u_from_rho_s] only implemented for eos 10') end select diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index 3632de54f..eccacacfa 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -136,6 +136,37 @@ subroutine get_eos_pressure_temp_mesa(den,eint,pres,temp) end subroutine get_eos_pressure_temp_mesa +!---------------------------------------------------------------- NOT TESTED YET! +!+ +! 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 + +!---------------------------------------------------------------- NOT TESTED YET! +!+ +! 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 c3ba9fbbd..25884fd0a 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -63,6 +63,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 @@ -552,7 +554,6 @@ end subroutine read_eos_mesa_gr ! 9. logS 10. dlnT/dlnP|S 11. Gamma1 12. gamma ! Note: ivout=1,2,3,4 returns the unlogged quantity -!!! I don't think a GR version of this is needed, the interpolation is the same. pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) real, intent(in) :: rho, eint real, intent(out) :: vout @@ -599,6 +600,56 @@ 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 + +pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) + 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) + logrho = log10(rho) + + ! Get the s and rho indices for looking up the tables + ns = 1 + int((logs - mesa_eos_s1) / mesa_eos_ds) + nrho = 1 + int((logrho - mesa_eos_rho1) / mesa_eos_drho) + + ! Allow extrapolation + if (ns < 1) ns=1 + if (ns > mesa_eos_ns-1) ns=mesa_eos_ns-1 + if (nrho < 1) nrho=1 + if (nrho > mesa_eos_nrho-1) nrho=mesa_eos_nrho-1 + + ds = (logs - mesa_eos_logSs(ns)) / mesa_eos_ds + drho = (logrho - mesa_eos_logRhos(nrho)) / mesa_eos_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 (I think? I can't tell if this is the correct way to do this) + if (ns > 1 .and. nrho > 1 .and. ns < mesa_eos_ns-1 .and. nrho < mesa_eos_nrho-1) then + call eos_cubic_spline_mesa_gr(ns,nrho,logs,logrho,ivout,vout,nx,dx) + if (ivout < 5) vout=10.d0**vout + if (present(ierr)) ierr = 0 + else + vout = 10.d0**((1.d0-ds) * (1.d0-drho) * mesa_ds_data(ns,nrho,ivout) + & + ds * (1.d0-drho) * mesa_ds_data(ns+1,nrho,ivout) + & + (1.d0-ds) * drho * mesa_ds_data(ns,nrho+1,ivout) + & + ds * drho * mesa_ds_data(ns+1,nrho+1,ivout)) + 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) @@ -655,6 +706,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_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_logss, mesa_eos_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_logss(s1-1) + x1=mesa_eos_logss(s1+0) + x2=mesa_eos_logss(s1+1) + x3=mesa_eos_logss(s1+2) + + tt=(s-x1)/mesa_eos_ds + + do j=1,4 + + y0 = mesa_ds_data(s1-1,rho1+j-2,n_var) + y1 = mesa_ds_data(s1+0,rho1+j-2,n_var) + y2 = mesa_ds_data(s1+1,rho1+j-2,n_var) + y3 = mesa_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_logrhos(rho1-1) + x1=mesa_eos_logrhos(rho1+0) + x2=mesa_eos_logrhos(rho1+1) + x3=mesa_eos_logrhos(rho1+2) + + y0=yy(1) + y1=yy(2) + y2=yy(3) + y3=yy(4) + + tt=(rho-x1)/mesa_eos_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 From 349b0d20aac8eef4df9f1df1150dafb546c186de Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 14:50:34 +1000 Subject: [PATCH 04/42] just saving for a second to compare with main branch --- src/main/eos.f90 | 12 ++++++++---- src/main/eos_mesa.f90 | 4 ++-- src/main/eos_mesa_microphysics.f90 | 5 +++-- 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index ce08d84d3..2910bb584 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -68,7 +68,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 @@ -1177,7 +1177,12 @@ end subroutine get_rho_from_p_s 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 +<<<<<<< HEAD use eos_idealplusrad, only:get_idealplusrad_tempfromrhoS +======= + use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_pres + use eos_mesa, only: get_eos_ptemp_from_rhos_mesa_gr +>>>>>>> d30372cf7 (just saving for a second to compare with main branch) use units, only:unit_density,unit_pressure,unit_ergg real, intent(in) :: S,mu,rho real, intent(inout) :: temp @@ -1223,10 +1228,11 @@ end subroutine get_p_from_rho_s subroutine get_u_from_rho_s(ieos,S,rho,u) use io, only:fatal use units, only:unit_density,unit_pressure,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,cgsp,cgss + real :: cgsrho,cgsp,cgss, cgsu ! change to cgs unit cgsrho = rho*unit_density @@ -1248,8 +1254,6 @@ subroutine get_u_from_rho_s(ieos,S,rho,u) end subroutine get_u_from_rho_s - - !----------------------------------------------------------------------- !+ ! Calculate mean molecular weight from X and Z, assuming complete diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index eccacacfa..ae5084594 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -61,8 +61,8 @@ subroutine init_eos_mesa(x,z,ierr) call read_eos_mesa(x,z,ierr) call get_opacity_constants_mesa call read_opacity_mesa(x,z) - mesa_eos_prefix="output_rhos_" - call read_eos_gr_mesa(x,z,ierr) ! read GR tables + mesa_eos_gr_prefix="output_rhos_" + call read_eos_mesa_gr(x,z,ierr) ! read GR tables end subroutine init_eos_mesa diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 25884fd0a..91710e94c 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -369,8 +369,9 @@ subroutine read_eos_mesa(x,z,ierr) ! Find the EoS tables filename = trim(mesa_eos_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' + print *, 'RAW filename: ', trim(filename) filename = find_phantom_datafile(filename,'eos/mesa') - + print *, 'RESOLVED filename: ', trim(filename) ! Note that the data exists mesa_eos_data_exists(i,j)=1 @@ -457,7 +458,7 @@ subroutine read_eos_mesa_gr(x,z,ierr) fmt1 = '(F4.2)' - fnum = 120 + fnum = 124 ierr = 0 ! following lines to prevent compiler warnings From 4286144b1c75398cbf49846fbbcaeccd0b99840a Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 15:22:33 +1000 Subject: [PATCH 05/42] cons2prim moved because I didn't commit it, so I had to rewrite it, sadly I don't have fitz to double check it rn, but I think it should work. Also, I still need to figure out why my tables aren't being read. --- src/main/cons2primsolver.f90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 7a616d1f7..b24d49117 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -198,6 +198,9 @@ 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,ent,dens,u) case (12) cgsdens = dens * unit_density cgsu = 1.5*rg*temp/gmw + radconst*temp**4/cgsdens @@ -261,6 +264,9 @@ 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,ent,dens,u) case (12) cgsdens = dens * unit_density cgsu = 1.5*rg*temp/gmw + radconst*temp**4/cgsdens From bf525e4e6a79c6ec9a85b44bcf8d317c8c6284b6 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 15:25:11 +1000 Subject: [PATCH 06/42] small bug --- src/main/cons2primsolver.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index b24d49117..98e413ff5 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -200,7 +200,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie select case(ieos) case (10) ! inputs and outputs are all in code units - call get_u_from_rho_s(ieos,ent,dens,u) + 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 @@ -266,7 +266,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie select case(ieos) case (10) ! inputs and outputs are all in code units - call get_u_from_rho_s(ieos,ent,dens,u) + 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 From c17180973ee0d1f2f3ac7abbda1bdaac7d0b50ac Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 15:25:49 +1000 Subject: [PATCH 07/42] small bug --- src/main/cons2primsolver.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 98e413ff5..81f25b4da 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 From 33b0f73d2c4ff1e179e933f179495d88451eb3c2 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 15:43:50 +1000 Subject: [PATCH 08/42] added debug outputs. something wrong with fnum reading --- src/main/cons2primsolver.f90 | 2 +- src/main/eos_mesa_microphysics.f90 | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 81f25b4da..5c7ad5cec 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -215,7 +215,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 diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 91710e94c..a31044650 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -369,15 +369,17 @@ subroutine read_eos_mesa(x,z,ierr) ! Find the EoS tables filename = trim(mesa_eos_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' - print *, 'RAW filename: ', trim(filename) + write(*,*) 'RAW filename: ', trim(filename) filename = find_phantom_datafile(filename,'eos/mesa') - print *, 'RESOLVED filename: ', trim(filename) + write(*,*) 'RESOLVED filename: ', trim(filename) ! 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 + write(*,*) "OPENING UNIT ", fnum, " FILE ", trim(filename) + 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) @@ -483,6 +485,7 @@ subroutine read_eos_mesa_gr(x,z,ierr) ! 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 + write(*,*) "OPENING UNIT ", fnum, " FILE ", trim(filename) open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted') read(fnum) mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 read(fnum)(mesa_eos_logrhos(k),k=1,mesa_eos_nrho) From 911d256c483aa13ada6fad502eb0ff101fe7bc1e Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 16:21:44 +1000 Subject: [PATCH 09/42] added specific labels for every gr related variable, but reading tables still crashes. idk why yet --- src/main/eos_mesa_microphysics.f90 | 148 +++++++++++++++-------------- 1 file changed, 76 insertions(+), 72 deletions(-) diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index a31044650..d509d42bf 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -43,14 +43,16 @@ module mesa_microphysics real, allocatable :: mesa_de_data0(:,:,:,:,:) ! EOS tables in entropy and density for GR - integer :: mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 - real :: mesa_eos_rho1, mesa_eos_s1, mesa_eos_ds, mesa_eos_drho + 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_logss(:), mesa_eos_logrhos(:) + 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_ds_data(:,:,:) -! real, allocatable :: mesa_eos0(:,:,:,:) - real, allocatable :: mesa_ds_data0(:,:,:,:,:) + 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 @@ -325,13 +327,14 @@ subroutine get_eos_constants_mesa(ierr) allocate(mesa_de_data0(mesa_eos_nz,mesa_eos_nh,mesa_eos_ne,mesa_eos_nv,mesa_eos_nvar2)) ! allocate GR tables (again, I wonder if it's best to make another subroutine for this) if (ierr /= 0) return - read(fnum) mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 + read(fnum) mesa_eos_gr_ns, mesa_eos_gr_nrho, mesa_eos_gr_nvar2 close(fnum) - allocate(mesa_eos_logss(mesa_eos_ns), mesa_eos_logrhos(mesa_eos_nrho)) - allocate(mesa_eos_gr_data_exists(mesa_eos_nz,mesa_eos_nh)) - allocate(mesa_ds_data(mesa_eos_ns,mesa_eos_nrho,mesa_eos_gr_nvar2)) - allocate(mesa_ds_data0(mesa_eos_nz,mesa_eos_nh,mesa_eos_ns,mesa_eos_nrho,mesa_eos_gr_nvar2)) + 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 @@ -460,7 +463,7 @@ subroutine read_eos_mesa_gr(x,z,ierr) fmt1 = '(F4.2)' - fnum = 124 + fnum = 126 ierr = 0 ! following lines to prevent compiler warnings @@ -469,11 +472,11 @@ subroutine read_eos_mesa_gr(x,z,ierr) dx = 0. !Loop over files - do i=1,mesa_eos_nz - write (zz,fmt1) mesa_eos_z(i) ! Z value as string + do i=1,mesa_eos_gr_nz + write (zz,fmt1) mesa_eos_gr_z(i) ! Z value as string - do j=1,mesa_eos_nh - write (hh,fmt1) mesa_eos_h(j) ! X value as string + do j=1,mesa_eos_gr_nh + write (hh,fmt1) mesa_eos_gr_h(j) ! X value as string !!!!!!!!!! have to rename these variables to not conflict with the above subroutine, but I can't be bothered to do that rn, so just ignore the fact that they have the same names as the above subroutine !!!!!!!!!! ! Find the EoS tables filename = trim(mesa_eos_gr_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' @@ -487,12 +490,12 @@ subroutine read_eos_mesa_gr(x,z,ierr) ! k, l and m hold the values of V, Eint and the data respectively write(*,*) "OPENING UNIT ", fnum, " FILE ", trim(filename) open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted') - read(fnum) mesa_eos_ns, mesa_eos_nrho, mesa_eos_gr_nvar2 - read(fnum)(mesa_eos_logrhos(k),k=1,mesa_eos_nrho) - read(fnum)(mesa_eos_logss(l),l=1,mesa_eos_ns) - do k=1,mesa_eos_nrho - do l=1,mesa_eos_ns - read(fnum) (mesa_ds_data0(i,j,l,k,m),m=1,mesa_eos_gr_nvar2) + 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 @@ -501,9 +504,9 @@ subroutine read_eos_mesa_gr(x,z,ierr) enddo ! Chooses the Z values from the available values which is greater than the input Z - nz2=mesa_eos_nz - do i=2,mesa_eos_nz-1 - if (mesa_eos_z(i) >= z) then + 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 @@ -511,41 +514,41 @@ subroutine read_eos_mesa_gr(x,z,ierr) ! 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_z(nz1))/(mesa_eos_z(nz2)-mesa_eos_z(nz1)) + 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_nh - if (mesa_eos_h(j) > x.or.j==mesa_eos_nh) then + 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_nh.and.mesa_eos_gr_data_exists(i,mesa_eos_nh-1)==0) then - nx2=mesa_eos_nh-1 + 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_h(nx1))/(mesa_eos_h(nx2)-mesa_eos_h(nx1)) + 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_ns - do k=1,mesa_eos_nrho + do l=1,mesa_eos_gr_ns + do k=1,mesa_eos_gr_nrho do m=1,mesa_eos_gr_nvar2 - mesa_ds_data(l,k,m)=(1.d0-dx)*(1.d0-dz)*mesa_ds_data0(nz1,nx1,l,k,m)+dz*(1.d0-dx)*mesa_ds_data0(nz2,nx1,l,k,m) & - +dx*(1.d0-dz)*mesa_ds_data0(nz1,nx2,l,k,m)+dx*dz*mesa_ds_data0(nz2,nx2,l,k,m) + 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_drho=mesa_eos_logrhos(2)-mesa_eos_logrhos(1) - mesa_eos_ds=mesa_eos_logss(2)-mesa_eos_logss(1) - mesa_eos_rho1=mesa_eos_logrhos(1) - mesa_eos_s1=mesa_eos_logss(1) + 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_eos0) - !deallocate(mesa_ds_data0) + !deallocate(mesa_eos_gr0) + !deallocate(mesa_gr_ds_data0) return @@ -623,29 +626,29 @@ pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) logrho = log10(rho) ! Get the s and rho indices for looking up the tables - ns = 1 + int((logs - mesa_eos_s1) / mesa_eos_ds) - nrho = 1 + int((logrho - mesa_eos_rho1) / mesa_eos_drho) + 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_ns-1) ns=mesa_eos_ns-1 + if (ns > mesa_eos_gr_ns-1) ns=mesa_eos_gr_ns-1 if (nrho < 1) nrho=1 - if (nrho > mesa_eos_nrho-1) nrho=mesa_eos_nrho-1 + if (nrho > mesa_eos_gr_nrho-1) nrho=mesa_eos_gr_nrho-1 - ds = (logs - mesa_eos_logSs(ns)) / mesa_eos_ds - drho = (logrho - mesa_eos_logRhos(nrho)) / mesa_eos_drho + 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 (I think? I can't tell if this is the correct way to do this) - if (ns > 1 .and. nrho > 1 .and. ns < mesa_eos_ns-1 .and. nrho < mesa_eos_nrho-1) then + 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 < 5) vout=10.d0**vout if (present(ierr)) ierr = 0 else - vout = 10.d0**((1.d0-ds) * (1.d0-drho) * mesa_ds_data(ns,nrho,ivout) + & - ds * (1.d0-drho) * mesa_ds_data(ns+1,nrho,ivout) + & - (1.d0-ds) * drho * mesa_ds_data(ns,nrho+1,ivout) + & - ds * drho * mesa_ds_data(ns+1,nrho+1,ivout)) + vout = 10.d0**((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 (present(ierr)) ierr = 1 ! warn if extrapolating endif @@ -711,12 +714,12 @@ 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_ds_data instead of mesa_de_data +! 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_logss, mesa_eos_logrhos +! use mesa_eos_gr_logss, mesa_eos_gr_logrhos implicit none integer, intent(in) :: s1, rho1, h1, n_var @@ -730,36 +733,36 @@ pure subroutine eos_cubic_spline_mesa_gr(s1,rho1,s,rho,n_var,z,h1,dh) real :: yy(4) integer :: j - x0=mesa_eos_logss(s1-1) - x1=mesa_eos_logss(s1+0) - x2=mesa_eos_logss(s1+1) - x3=mesa_eos_logss(s1+2) + 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_ds + tt=(s-x1)/mesa_eos_gr_ds do j=1,4 - y0 = mesa_ds_data(s1-1,rho1+j-2,n_var) - y1 = mesa_ds_data(s1+0,rho1+j-2,n_var) - y2 = mesa_ds_data(s1+1,rho1+j-2,n_var) - y3 = mesa_ds_data(s1+2,rho1+j-2,n_var) + 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_logrhos(rho1-1) - x1=mesa_eos_logrhos(rho1+0) - x2=mesa_eos_logrhos(rho1+1) - x3=mesa_eos_logrhos(rho1+2) + 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_drho + 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) @@ -805,11 +808,12 @@ subroutine deallocate_arrays_mesa if (allocated(mesa_de_data)) deallocate(mesa_de_data) if (allocated(mesa_de_data0)) deallocate(mesa_de_data0) !eos gr (I wonder if it's best to make another subroutine for deallocating the GR tables) - if (allocated(mesa_eos_logss)) deallocate(mesa_eos_logss) - if (allocated(mesa_eos_logrhos)) deallocate(mesa_eos_logrhos) + 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_ds_data)) deallocate(mesa_ds_data) - if (allocated(mesa_ds_data0)) deallocate(mesa_ds_data0) + 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 From cd73c1f07643b2fc23e11403e00cd20e687820ca Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 16:57:03 +1000 Subject: [PATCH 10/42] added get_eos_constants_mesa_gr because the crashes won't stop. sitll not tested --- src/main/eos_mesa_microphysics.f90 | 54 ++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index d509d42bf..488833f86 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -340,6 +340,60 @@ subroutine get_eos_constants_mesa(ierr) 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 (again, I wonder if it's best to make another subroutine for this) + 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 From 89a7186e88dab6355443ceb36c309f6ccdbdbdcc Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 17:00:28 +1000 Subject: [PATCH 11/42] minor bug --- src/main/eos_mesa.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index ae5084594..9aa3b19d1 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -57,7 +57,8 @@ subroutine init_eos_mesa(x,z,ierr) call get_eos_constants_mesa(ierr) if (ierr /= 0) return - + call get_eos_constants_mesa_gr(ierr) + if (ierr /= 0) return call read_eos_mesa(x,z,ierr) call get_opacity_constants_mesa call read_opacity_mesa(x,z) From a651bcfc7de4d8bfba4c416497d7cd4a43764c77 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 17:32:32 +1000 Subject: [PATCH 12/42] tables are being read now! there is a new bug which I will fix later about GR only working for ieos2,12, or 11 --- src/main/eos_mesa.f90 | 3 ++- src/main/eos_mesa_microphysics.f90 | 10 ---------- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index 9aa3b19d1..e0604f000 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -57,12 +57,13 @@ subroutine init_eos_mesa(x,z,ierr) call get_eos_constants_mesa(ierr) if (ierr /= 0) return + mesa_eos_gr_prefix="output_rhos_" + call get_eos_constants_mesa_gr(ierr) if (ierr /= 0) return call read_eos_mesa(x,z,ierr) call get_opacity_constants_mesa call read_opacity_mesa(x,z) - mesa_eos_gr_prefix="output_rhos_" call read_eos_mesa_gr(x,z,ierr) ! read GR tables end subroutine init_eos_mesa diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 488833f86..e9c40af5d 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -325,16 +325,6 @@ subroutine get_eos_constants_mesa(ierr) allocate(mesa_de_data(mesa_eos_ne,mesa_eos_nv,mesa_eos_nvar2)) 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)) -! allocate GR tables (again, I wonder if it's best to make another subroutine for this) - 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 From f5a0a1041dc0a0b17aab8b9fe971ac0b8ffe0653 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Sat, 13 Jun 2026 23:23:43 +1000 Subject: [PATCH 13/42] added eos=10 to eos.f90 cases for gr, also removed some comments. now there is an issue probably with the tables themselves which I'll need to test --- src/main/eos.f90 | 2 +- src/main/eos_mesa_microphysics.f90 | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 2910bb584..e9499f62c 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -164,7 +164,7 @@ 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',& var='eos_type',val=real(eos_type)) diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index e9c40af5d..644f7aa88 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -416,16 +416,13 @@ subroutine read_eos_mesa(x,z,ierr) ! Find the EoS tables filename = trim(mesa_eos_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' - write(*,*) 'RAW filename: ', trim(filename) filename = find_phantom_datafile(filename,'eos/mesa') - write(*,*) 'RESOLVED filename: ', trim(filename) ! 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 - write(*,*) "OPENING UNIT ", fnum, " FILE ", trim(filename) open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted') read(fnum) mesa_eos_ne, mesa_eos_nv, mesa_eos_nvar2 @@ -532,7 +529,6 @@ subroutine read_eos_mesa_gr(x,z,ierr) ! 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 - write(*,*) "OPENING UNIT ", fnum, " FILE ", trim(filename) 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) From d78d5ae23a1ded4d91dd43a1dcbb0616838bf013 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 1 Jul 2026 17:16:55 +1000 Subject: [PATCH 14/42] minor addition --- src/main/eos.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index e9499f62c..eb0324e26 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1116,6 +1116,7 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in) end function entropy +! input is in cgs units, output is in cgss/(kb*unit_ergg) real function get_entropy(rho,pres,mu_in,ieos) use units, only:unit_density,unit_pressure,unit_ergg use physcon, only:kboltz From e735f3d430e08def16fafe63a4600660059f7c46 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 1 Jul 2026 18:39:54 +1000 Subject: [PATCH 15/42] corrected the units for entropy when going back and forth from the mesa tables, because they are specific entropy divided by avo*kboltz --- src/main/eos.f90 | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index eb0324e26..62ba8884f 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -164,7 +164,7 @@ 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,10,11,12/)==eos_type)) then + if (gr .and. .not.any((/2,4,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',& var='eos_type',val=real(eos_type)) @@ -1116,7 +1116,6 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in) end function entropy -! input is in cgs units, output is in cgss/(kb*unit_ergg) real function get_entropy(rho,pres,mu_in,ieos) use units, only:unit_density,unit_pressure,unit_ergg use physcon, only:kboltz @@ -1178,12 +1177,8 @@ end subroutine get_rho_from_p_s 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 -<<<<<<< HEAD - use eos_idealplusrad, only:get_idealplusrad_tempfromrhoS -======= use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_pres use eos_mesa, only: get_eos_ptemp_from_rhos_mesa_gr ->>>>>>> d30372cf7 (just saving for a second to compare with main branch) use units, only:unit_density,unit_pressure,unit_ergg real, intent(in) :: S,mu,rho real, intent(inout) :: temp From 8a50944c64667f3b24a7d5d62e03bb7647e1344b Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 1 Jul 2026 18:41:36 +1000 Subject: [PATCH 16/42] minor debug --- src/main/eos.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 62ba8884f..f1574fb15 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1223,7 +1223,8 @@ end subroutine get_p_from_rho_s !----------------------------------------------------------------------- subroutine get_u_from_rho_s(ieos,S,rho,u) use io, only:fatal - use units, only:unit_density,unit_pressure,unit_ergg + use units, only:unit_density,unit_pressure,unit_ergg,umass + use physcon, only:kboltz,avogadro use eos_mesa, only: get_eos_u_from_rhos_mesa_gr real, intent(in) :: S,rho real, intent(out) :: u From f7bba6c9ca986d542e5cbdba5bd34e8c1d9b8bea Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 1 Jul 2026 19:12:10 +1000 Subject: [PATCH 17/42] I removed the k_b denom from entropy definition. eos mesa tables back and forth was getting insane --- src/main/eos.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index f1574fb15..25a029005 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1116,6 +1116,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 @@ -1133,8 +1134,8 @@ real function get_entropy(rho,pres,mu_in,ieos) case default cgss = entropy(cgsrho,cgspres,mu_in,1) end select - cgss = cgss/kboltz ! s/kb - get_entropy = cgss/unit_ergg +! cgss = cgss/kboltz ! s/kb + get_entropy = cgss/unit_ergg ! units in erg/grK, here it turns to code units end function get_entropy @@ -1189,8 +1190,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) ! change to cgs unit cgsrho = rho*unit_density - cgss = s*unit_ergg - if (present(niter_out)) niter_out = 0 + cgss = S*unit_ergg select case (ieos) case (2,5) @@ -1223,7 +1223,7 @@ end subroutine get_p_from_rho_s !----------------------------------------------------------------------- subroutine get_u_from_rho_s(ieos,S,rho,u) use io, only:fatal - use units, only:unit_density,unit_pressure,unit_ergg,umass + use units, only:unit_density,unit_pressure,unit_ergg use physcon, only:kboltz,avogadro use eos_mesa, only: get_eos_u_from_rhos_mesa_gr real, intent(in) :: S,rho @@ -1233,7 +1233,7 @@ subroutine get_u_from_rho_s(ieos,S,rho,u) ! change to cgs unit cgsrho = rho*unit_density - cgss = s*unit_ergg + cgss = S*unit_ergg select case (ieos) case(10) From 2db3074e1522b841b4e48257aeefb933e1eb5a5a Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 1 Jul 2026 20:25:34 +1000 Subject: [PATCH 18/42] made an important change . note that in eos mesa some tables are not in log, so you need to specify which are and areny --- src/main/eos.f90 | 10 ++++++---- src/main/eos_mesa_microphysics.f90 | 9 +++++---- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 25a029005..bd8675bc9 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1185,8 +1185,10 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) real, intent(inout) :: temp real, intent(out) :: P integer, intent(in) :: ieos - integer, intent(out), optional :: niter_out - real :: cgsrho,cgspres,cgss + real :: corr,df,f,temp_new,cgsrho,cgsp,cgss + real, parameter :: eoserr=1e-12 + integer :: niter + integer, parameter :: nitermax = 1000 ! change to cgs unit cgsrho = rho*unit_density @@ -1198,8 +1200,8 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) cgspres = cgsrho*Rg*temp / mu case(10) !!! For GR, not tested yet! + mesas = cgss/(kboltz*avogadro) ! the MESA tables are specific entropy divided by (avo*kerg). call get_eos_ptemp_from_rhos_mesa_gr(cgsrho,cgss,cgsP,temp) - case (12) call get_idealplusrad_tempfromrhoS(cgsrho,cgss,mu,temp,cgspres,niter_out) case default @@ -1238,7 +1240,7 @@ subroutine get_u_from_rho_s(ieos,S,rho,u) select case (ieos) case(10) !!! For GR, not tested yet! - call get_eos_u_from_rhos_mesa_gr(cgsrho,cgss,cgsu) + call get_eos_u_from_rhos_mesa_gr(cgsrho,cgss,cgsu) case default cgsu = 0. diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 644f7aa88..6546d42cb 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -633,7 +633,7 @@ pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) ! 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) 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) 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) + & @@ -651,8 +651,9 @@ end subroutine getvalue_mesa ! 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 @@ -662,7 +663,7 @@ pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) real :: dx integer :: nx - logs = log10(s) + 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 @@ -682,7 +683,7 @@ pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) ! 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) 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 < 5) vout=10.d0**vout + 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 = 10.d0**((1.d0-ds) * (1.d0-drho) * mesa_gr_ds_data(ns,nrho,ivout) + & From 7addc6d72fde413ed06b47f3c49bede3b775030f Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Thu, 2 Jul 2026 15:17:02 +1000 Subject: [PATCH 19/42] just adding a checkpoint here. the entropy is now being read and stored in eosvars but now has to somehow go to be used in relaxstar --- src/main/cons2primsolver.f90 | 13 ++++++++----- src/main/eos.f90 | 6 +++--- src/main/eos_mesa_microphysics.f90 | 2 +- src/main/part.F90 | 5 +++-- src/setup/relax_star.f90 | 1 + src/setup/set_star_utils.f90 | 19 ++++++++++++++++--- 6 files changed, 32 insertions(+), 14 deletions(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 5c7ad5cec..3caa66145 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -157,7 +157,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie logical :: converged,have_eos_cache real :: gcov(0:3,0:3) ierr = 0 - + write(*,*) 'ALI: DEBUG ieos=', ieos, ' ien_type=', ien_type ! Get metric components from metric array call unpack_metric(metrici,gcov=gcov,gammaijUP=gammaijUP,alpha=alpha,betadown=betadown,betaUP=betaUP) @@ -193,14 +193,15 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie p = (gamma-1.)*dens*polyk elseif (ien_type == ien_entropy_s) then call get_p_from_rho_s(ieos,en,dens,gmw,P,temp) - dens_eos = dens - p_eos = P - temp_eos = temp - have_eos_cache = .true. + ! write(*,*) 'ALI, en = ',en,' dens = ',dens,' P = ',P,' temp = ',temp + ! write(*,*) 'ENTER entropy branch' + ! write(*,*) 'ieos=', ieos select case(ieos) case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) + write(*,*) 'ALI 1, went into case 10 EOS MESA tables' + case (12) cgsdens = dens * unit_density cgsu = 1.5*rg*temp/gmw + radconst*temp**4/cgsdens @@ -267,6 +268,8 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) + write(*,*) 'ALI 2, went into case 10 EOS MESA tables' + case (12) cgsdens = dens * unit_density cgsu = 1.5*rg*temp/gmw + radconst*temp**4/cgsdens diff --git a/src/main/eos.f90 b/src/main/eos.f90 index bd8675bc9..e26d1bf07 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1123,7 +1123,8 @@ real function get_entropy(rho,pres,mu_in,ieos) integer, intent(in) :: ieos real, intent(in) :: rho,pres,mu_in real :: cgsrho,cgspres,cgss - + +! write(*,*) 'ALI: GET_ENTROPY Called: ieos = ',ieos cgsrho = rho * unit_density cgspres = pres * unit_pressure select case (ieos) @@ -1134,7 +1135,7 @@ real function get_entropy(rho,pres,mu_in,ieos) case default cgss = entropy(cgsrho,cgspres,mu_in,1) end select -! cgss = cgss/kboltz ! s/kb + cgss = cgss/kboltz ! s/kb get_entropy = cgss/unit_ergg ! units in erg/grK, here it turns to code units end function get_entropy @@ -1200,7 +1201,6 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) cgspres = cgsrho*Rg*temp / mu case(10) !!! For GR, not tested yet! - mesas = cgss/(kboltz*avogadro) ! the MESA tables are specific entropy divided by (avo*kerg). call get_eos_ptemp_from_rhos_mesa_gr(cgsrho,cgss,cgsP,temp) case (12) call get_idealplusrad_tempfromrhoS(cgsrho,cgss,mu,temp,cgspres,niter_out) diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 6546d42cb..ea2ba2c77 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -633,7 +633,7 @@ pure subroutine getvalue_mesa(rho,eint,ivout,vout,ierr) ! 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) 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. ! this is applied only for the first four columns of the data which are stored in log10 form. + if (ivout < 5) 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) + & diff --git a/src/main/part.F90 b/src/main/part.F90 index b7b333cc2..3ebbef188 100644 --- a/src/main/part.F90 +++ b/src/main/part.F90 @@ -159,9 +159,10 @@ module part iX = 5, & iZ = 6, & igamma = 7, & - maxeosvars = 7 + ientr = 8, & + maxeosvars = 8 character(len=*), parameter :: eos_vars_label(maxeosvars) = & - (/'pressure ','sound speed','temperature','mu ','H fraction ','metallicity','gamma '/) + (/'pressure ','sound speed','temperature','mu ','H fraction ','metallicity','gamma ','entropy '/) ! !--energy_variables diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index 97283e5f8..62de6dd0d 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -273,6 +273,7 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& print "(a,i0,a,i0,a,2pf6.2,2(a,1pg11.3))",& ' Relaxing star: Iter ',nits,'/',maxits, & ', dens error:',rmserr,'%, R*:',rmax,' Ekin/Epot:',ekin/abs(epot) + write(*,*) 'Ali: maxvxyzu = ', maxvxyzu, ' ieos = ', ieos, ' ieos_prev = ', ieos_prev endif endif ! diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 241622a2a..8138c1162 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -465,9 +465,10 @@ end subroutine set_star_composition !----------------------------------------------------------------------- subroutine set_star_thermalenergy(ieos,den,pres,temp,r,npts,npart,xyzh,vxyzu,rad,eos_vars,& relaxed,use_var_comp,initialtemp,polyk_in,npin,x0) - use part, only:rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi,icv,& - aprmassoftype,apr_level,radprop - use eos, only:equationofstate,calc_temp_and_ene,eos_outputs_mu,get_cv,gmw + use part, only:rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi,icv,ientrop,& + aprmassoftype,apr_level,radprop,ien_type,ien_entropy_s ! added some stuff ALI + use eos, only:equationofstate,calc_temp_and_ene,eos_outputs_mu,get_cv,gmw,& + get_entropy use radiation_utils, only:radxi_from_Trad use table_utils, only:yinterp use units, only:unit_density,unit_ergg,unit_pressure @@ -579,6 +580,18 @@ subroutine set_star_thermalenergy(ieos,den,pres,temp,r,npts,npart,xyzh,vxyzu,rad endif eos_vars(itemp,i) = tempi vxyzu(4,i) = eni / unit_ergg + + ! ALI: ADDED: set entropy for ieos=10 EOS MESA for GR, under construction + ! mu won't be used by entropy( for eos_mesa anyway + if (ien_type == ien_entropy_s) then + if (use_var_comp) then + mu = eos_vars(imu,i) + else + mu = gmw + endif + eos_vars(ientrop,i) = get_entropy(densi,presi,mu,ieos) + endif + end select endif enddo From 9f303c6e9424fe197c6d5772d296d9fd52fd49a5 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Fri, 3 Jul 2026 15:12:24 +1000 Subject: [PATCH 20/42] calling entropy and thermal energy in relax_star to be from mesa tables for eosmesa, this well affect classical star relaxing as well as grstar and make them dependent on entropy --- src/main/eos.f90 | 35 ++++++++++++++++++++++++++++++++++- src/setup/relax_star.f90 | 30 ++++++++++++++++++++++-------- 2 files changed, 56 insertions(+), 9 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index e26d1bf07..745d78089 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -69,7 +69,7 @@ module eos 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,get_u_from_rho_s - public :: calc_rho_from_PT,get_entropy,get_p_from_rho_s + public :: calc_rho_from_PT,get_entropy,get_p_from_rho_s,get_entropy_vec, get_u_from_rho_s_vec public :: init_eos,finish_eos public :: write_options_eos,read_options_eos,set_defaults_eos public :: write_headeropts_eos,read_headeropts_eos @@ -1140,6 +1140,27 @@ real function get_entropy(rho,pres,mu_in,ieos) end function get_entropy +! vectorised version of get_entropy for use in relax_star +subroutine get_entropy_vec(rho, pres, mu_in, ieos, S) + use io, only: fatal + integer, intent(in) :: ieos + real, intent(in) :: rho(:), pres(:) + real, intent(in) :: mu_in + real, intent(out) :: S(:) + + integer :: i, n + + n = size(rho) + + if (size(pres) /= n .or. size(S) /= n) then + call fatal('get_entropy_vec','array size mismatch') + endif + + do i = 1, n + S(i) = get_entropy(rho(i), pres(i), mu_in, ieos) + end do + +end subroutine get_entropy_vec !----------------------------------------------------------------------- !+ ! Calculate density given pressure and entropy using Newton-Raphson @@ -1252,6 +1273,18 @@ subroutine get_u_from_rho_s(ieos,S,rho,u) end subroutine get_u_from_rho_s +! vectorised version of get_u_from_rho_s for use in relax_star +subroutine get_u_from_rho_s_vec(ieos,S,rho,u,n) + integer, intent(in) :: ieos, n + real, intent(in) :: S(n), rho(n) + real, intent(out) :: u(n) + + integer :: i + + do i = 1, n + call get_u_from_rho_s(ieos, S(i), rho(i), u(i)) + end do +end subroutine get_u_from_rho_s_vec !----------------------------------------------------------------------- !+ diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index 62de6dd0d..af561b20b 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -69,7 +69,8 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& use io, only:error,warning,fatal,id,master use fileutils, only:getnextfilename use readwrite_dumps, only:write_fulldump - use eos, only:gamma,eos_outputs_mu,ieos,eos_has_pressure_without_u,polyk + use eos, only:gamma,eos_outputs_mu,ieos,eos_has_pressure_without_u,polyk,& + get_entropy_vec, get_u_from_rho_s_vec !Ali use physcon, only:pi use options, only:iexternalforce use io_summary, only:summary_initialise @@ -95,6 +96,8 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& logical, parameter :: fix_entrop = .true. ! fix entropy instead of thermal energy logical :: write_files character(len=20) :: filename,mylabel + real :: gmw_dummy ! Ali + gmw_dummy = 1.0 ! Ali i1 = 0 if (present(npin)) i1 = npin ! starting position in particle array @@ -176,13 +179,24 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& ! define utherm(r) based on P(r) and rho(r) ! and use this to set the thermal energy of all particles ! - where (rho > epsilon(0.) .and. gamma > 1.) - entrop = pr/rho**gamma - utherm = pr/(rho*(gamma-1.)) - elsewhere - entrop = 0. - utherm = 0. - end where + write(*,*) 'Ali before entropy set: ieos = ', ieos + if (ieos == 10) then + + call get_entropy_vec(rho, pr, gmw_dummy, ieos, entrop) + call get_u_from_rho_s_vec(ieos, entrop, rho, utherm, nt) + + else + + where (rho > epsilon(0.) .and. gamma > 1.) + entrop = pr / rho**gamma + utherm = pr / (rho*(gamma-1.)) + elsewhere + entrop = 0. + utherm = 0. + end where + + endif + if (any(utherm(1:nt-1) <= 0.) .and. .not.eos_has_pressure_without_u(ieos)) then call error('relax_star','relax-o-matic needs non-zero pressure array set in order to work') From c0700b69d29c380f4a17f7a1d9d2c3dcf4051438 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Fri, 3 Jul 2026 16:55:36 +1000 Subject: [PATCH 21/42] added eos mesa to relax star, now there will be a lot of debugging, also to see whehter Daniel is happy with this at all --- src/main/cons2primsolver.f90 | 6 +++--- src/setup/relax_star.f90 | 13 ++++++++----- src/setup/set_star_utils.f90 | 19 +++---------------- 3 files changed, 14 insertions(+), 24 deletions(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 3caa66145..f89de0f46 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -157,7 +157,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie logical :: converged,have_eos_cache real :: gcov(0:3,0:3) ierr = 0 - write(*,*) 'ALI: DEBUG ieos=', ieos, ' ien_type=', ien_type +! write(*,*) 'ALI: DEBUG ieos=', ieos, ' ien_type=', ien_type ! Get metric components from metric array call unpack_metric(metrici,gcov=gcov,gammaijUP=gammaijUP,alpha=alpha,betadown=betadown,betaUP=betaUP) @@ -200,7 +200,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) - write(*,*) 'ALI 1, went into case 10 EOS MESA tables' + ! write(*,*) 'ALI 1, went into case 10 EOS MESA tables' case (12) cgsdens = dens * unit_density @@ -268,7 +268,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) - write(*,*) 'ALI 2, went into case 10 EOS MESA tables' + ! write(*,*) 'ALI 2, went into case 10 EOS MESA tables' case (12) cgsdens = dens * unit_density diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index af561b20b..cbf6e834b 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -179,11 +179,11 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& ! define utherm(r) based on P(r) and rho(r) ! and use this to set the thermal energy of all particles ! - write(*,*) 'Ali before entropy set: ieos = ', ieos - if (ieos == 10) then + write(*,*) 'Ali before entropy set: ieos_prev = ', ieos_prev + if (ieos_prev == 10) then - call get_entropy_vec(rho, pr, gmw_dummy, ieos, entrop) - call get_u_from_rho_s_vec(ieos, entrop, rho, utherm, nt) + call get_entropy_vec(rho, pr, gmw_dummy, ieos_prev, entrop) + call get_u_from_rho_s_vec(ieos_prev, entrop, rho, utherm, nt) else @@ -545,7 +545,10 @@ subroutine set_options_for_relaxation(tdyn) ! ! turn on settings appropriate to relaxation ! - if (maxvxyzu >= 4 .and. .not.eos_has_pressure_without_u(ieos)) ieos = 2 + if (ieos /= 10) then + if (maxvxyzu >= 4 .and. .not.eos_has_pressure_without_u(ieos)) ieos = 2 + endif +! if (maxvxyzu >= 4 .and. .not.eos_has_pressure_without_u(ieos)) ieos = 2 if (tdyn > 0.) then idamp = 2 tdyn_s = tdyn*utime diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 8138c1162..241622a2a 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -465,10 +465,9 @@ end subroutine set_star_composition !----------------------------------------------------------------------- subroutine set_star_thermalenergy(ieos,den,pres,temp,r,npts,npart,xyzh,vxyzu,rad,eos_vars,& relaxed,use_var_comp,initialtemp,polyk_in,npin,x0) - use part, only:rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi,icv,ientrop,& - aprmassoftype,apr_level,radprop,ien_type,ien_entropy_s ! added some stuff ALI - use eos, only:equationofstate,calc_temp_and_ene,eos_outputs_mu,get_cv,gmw,& - get_entropy + use part, only:rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi,icv,& + aprmassoftype,apr_level,radprop + use eos, only:equationofstate,calc_temp_and_ene,eos_outputs_mu,get_cv,gmw use radiation_utils, only:radxi_from_Trad use table_utils, only:yinterp use units, only:unit_density,unit_ergg,unit_pressure @@ -580,18 +579,6 @@ subroutine set_star_thermalenergy(ieos,den,pres,temp,r,npts,npart,xyzh,vxyzu,rad endif eos_vars(itemp,i) = tempi vxyzu(4,i) = eni / unit_ergg - - ! ALI: ADDED: set entropy for ieos=10 EOS MESA for GR, under construction - ! mu won't be used by entropy( for eos_mesa anyway - if (ien_type == ien_entropy_s) then - if (use_var_comp) then - mu = eos_vars(imu,i) - else - mu = gmw - endif - eos_vars(ientrop,i) = get_entropy(densi,presi,mu,ieos) - endif - end select endif enddo From 43e573cc95b27c1b2cf33e8b21370bc37e4c34b2 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Mon, 6 Jul 2026 08:48:16 +1000 Subject: [PATCH 22/42] upon consulting ryo, we decided not to add entropy of mesa during relaxation because the entropy at that stage shouldn't be tampered with --- src/main/part.F90 | 5 ++--- src/setup/relax_star.f90 | 36 +++++++++--------------------------- 2 files changed, 11 insertions(+), 30 deletions(-) diff --git a/src/main/part.F90 b/src/main/part.F90 index 3ebbef188..b7b333cc2 100644 --- a/src/main/part.F90 +++ b/src/main/part.F90 @@ -159,10 +159,9 @@ module part iX = 5, & iZ = 6, & igamma = 7, & - ientr = 8, & - maxeosvars = 8 + maxeosvars = 7 character(len=*), parameter :: eos_vars_label(maxeosvars) = & - (/'pressure ','sound speed','temperature','mu ','H fraction ','metallicity','gamma ','entropy '/) + (/'pressure ','sound speed','temperature','mu ','H fraction ','metallicity','gamma '/) ! !--energy_variables diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index cbf6e834b..97283e5f8 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -69,8 +69,7 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& use io, only:error,warning,fatal,id,master use fileutils, only:getnextfilename use readwrite_dumps, only:write_fulldump - use eos, only:gamma,eos_outputs_mu,ieos,eos_has_pressure_without_u,polyk,& - get_entropy_vec, get_u_from_rho_s_vec !Ali + use eos, only:gamma,eos_outputs_mu,ieos,eos_has_pressure_without_u,polyk use physcon, only:pi use options, only:iexternalforce use io_summary, only:summary_initialise @@ -96,8 +95,6 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& logical, parameter :: fix_entrop = .true. ! fix entropy instead of thermal energy logical :: write_files character(len=20) :: filename,mylabel - real :: gmw_dummy ! Ali - gmw_dummy = 1.0 ! Ali i1 = 0 if (present(npin)) i1 = npin ! starting position in particle array @@ -179,24 +176,13 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& ! define utherm(r) based on P(r) and rho(r) ! and use this to set the thermal energy of all particles ! - write(*,*) 'Ali before entropy set: ieos_prev = ', ieos_prev - if (ieos_prev == 10) then - - call get_entropy_vec(rho, pr, gmw_dummy, ieos_prev, entrop) - call get_u_from_rho_s_vec(ieos_prev, entrop, rho, utherm, nt) - - else - - where (rho > epsilon(0.) .and. gamma > 1.) - entrop = pr / rho**gamma - utherm = pr / (rho*(gamma-1.)) - elsewhere - entrop = 0. - utherm = 0. - end where - - endif - + where (rho > epsilon(0.) .and. gamma > 1.) + entrop = pr/rho**gamma + utherm = pr/(rho*(gamma-1.)) + elsewhere + entrop = 0. + utherm = 0. + end where if (any(utherm(1:nt-1) <= 0.) .and. .not.eos_has_pressure_without_u(ieos)) then call error('relax_star','relax-o-matic needs non-zero pressure array set in order to work') @@ -287,7 +273,6 @@ subroutine relax_star(nt,rho,pr,temp,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,& print "(a,i0,a,i0,a,2pf6.2,2(a,1pg11.3))",& ' Relaxing star: Iter ',nits,'/',maxits, & ', dens error:',rmserr,'%, R*:',rmax,' Ekin/Epot:',ekin/abs(epot) - write(*,*) 'Ali: maxvxyzu = ', maxvxyzu, ' ieos = ', ieos, ' ieos_prev = ', ieos_prev endif endif ! @@ -545,10 +530,7 @@ subroutine set_options_for_relaxation(tdyn) ! ! turn on settings appropriate to relaxation ! - if (ieos /= 10) then - if (maxvxyzu >= 4 .and. .not.eos_has_pressure_without_u(ieos)) ieos = 2 - endif -! if (maxvxyzu >= 4 .and. .not.eos_has_pressure_without_u(ieos)) ieos = 2 + if (maxvxyzu >= 4 .and. .not.eos_has_pressure_without_u(ieos)) ieos = 2 if (tdyn > 0.) then idamp = 2 tdyn_s = tdyn*utime From 78ed890939a4d8dc7b3765e7139decf3bbc59d07 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Mon, 6 Jul 2026 19:04:17 +1000 Subject: [PATCH 23/42] some useful debugging flags I want to save --- src/main/cons2primsolver.f90 | 28 +++++++++++++++++++++++++--- src/main/eos.f90 | 27 ++++++++++++++++++++++----- 2 files changed, 47 insertions(+), 8 deletions(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index f89de0f46..1719a5477 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -78,6 +78,7 @@ subroutine primitive2conservative(x,metrici,v,dens,u,P,rho,pmom,en,ien_type) use metric_tools, only:unpack_metric use io, only:error use eos, only:gmw,get_entropy + use units, only:unit_density,unit_ergg ! ali added real, intent(in) :: x(1:3),metrici(:,:,:) real, intent(in) :: dens,v(1:3),u,P real, intent(out) :: rho,pmom(1:3),en @@ -113,7 +114,12 @@ subroutine primitive2conservative(x,metrici,v,dens,u,P,rho,pmom,en,ien_type) if (ien_type == ien_etotal) then en = U0*enth*gvv + (1.+u)/U0 elseif (ien_type == ien_entropy_s) then + ! write(*,'(A,ES15.7)') 'before get_entropy in cons2primsolver:entropy (code units): ', en en = get_entropy(dens,P,gmw,ieos) + ! write(*,'(A,ES15.7)') 'after get_entropy in cons2primsolver:entropy (code units): ', en + ! write(*,'(A,ES15.7)') 'get_entropy:rho (cgs): ', dens*unit_density + ! write(*,'(A,ES15.7)') 'get_entropy:P (cgs): ', P*unit_ergg/unit_density + ! write(*,'(A,ES15.7)') 'get_entropy:u (cgs): ', u*unit_ergg else if (u > 0) then gam1 = 1. + P/(dens*u) @@ -157,7 +163,13 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie logical :: converged,have_eos_cache real :: gcov(0:3,0:3) ierr = 0 -! write(*,*) 'ALI: DEBUG ieos=', ieos, ' ien_type=', ien_type +! write(*,*) 'L158: ALI cons2primsolver: DEBUG ieos=', ieos, ' ien_type=', ien_type +! write(*,'(A,ES15.7)') 'L158:rho (cgs): ', dens*unit_density +! write(*,'(A,ES15.7)') 'L158:entropy (cgs): ', en*unit_ergg +! write(*,'(A,ES15.7)') 'L158:u (cgs): ', u*unit_ergg +! write(*,'(A,ES15.7)') 'L158:P (cgs): ', P*unit_ergg/unit_density +! write(*,'(A,ES15.7)') 'L158:Temperature (K): ', temp + ! Get metric components from metric array call unpack_metric(metrici,gcov=gcov,gammaijUP=gammaijUP,alpha=alpha,betadown=betadown,betaUP=betaUP) @@ -200,7 +212,12 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) - ! write(*,*) 'ALI 1, went into case 10 EOS MESA tables' + ! write(*,*) 'L200: ALI cons2primsolver, went into case 10 EOS MESA tables' + ! write(*,'(A,ES15.7)') 'L200:rho (cgs): ', dens*unit_density + ! write(*,'(A,ES15.7)') 'L200:entropy (cgs): ', en*unit_ergg + ! write(*,'(A,ES15.7)') 'L200:u (cgs): ', u*unit_ergg + ! write(*,'(A,ES15.7)') 'L200:P (cgs): ', P*unit_ergg/unit_density + ! write(*,'(A,ES15.7)') 'L200:Temperature (K): ', temp case (12) cgsdens = dens * unit_density @@ -268,7 +285,12 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) - ! write(*,*) 'ALI 2, went into case 10 EOS MESA tables' + ! write(*,*) 'L263: ALI cons2primsolver, went into case 10 EOS MESA tables' + ! write(*,'(A,ES15.7)') 'L263:rho (cgs): ', dens*unit_density + ! write(*,'(A,ES15.7)') 'L263:entropy (cgs): ', en*unit_ergg + ! write(*,'(A,ES15.7)') 'L263:u (cgs): ', u*unit_ergg + ! write(*,'(A,ES15.7)') 'L263:P (cgs): ', P*unit_ergg/unit_density + ! write(*,'(A,ES15.7)') 'L263:Temperature (K): ', temp case (12) cgsdens = dens * unit_density diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 745d78089..b3d79c0eb 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -334,6 +334,18 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam if (present(gamma_local)) gamma_local = gam1 ! gamma is an output if (present(mu_local)) mu_local = 1./get_eos_1overmu_mesa(cgsrhoi,cgseni) if (ierr /= 0) call warning('eos_mesa','extrapolating off tables') + if (ierr /= 0) then + write(*,'(A)') '*** MESA EOS extrapolating off tables ***' + ! write(*,'(A,ES15.7)') 'rho (code units): ', rhoi + write(*,'(A,ES15.7)') 'rho (cgs): ', cgsrhoi + ! write(*,'(A,ES15.7)') 'u (code units): ', eni + write(*,'(A,ES15.7)') 'u (cgs): ', cgseni + ! write(*,'(A,ES15.7)') 'P (code units): ', presi + write(*,'(A,ES15.7)') 'P (cgs): ', cgspresi + write(*,'(A,F12.6)') 'gamma1: ', gam1 + write(*,'(A,ES15.7)') 'Temperature (K): ', temperaturei + call warning('eos_mesa','extrapolating off tables') + endif case(11) ! @@ -1101,14 +1113,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, + ! because for vars >5 in the tables the log becomes the actual value) 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)') @@ -1135,8 +1148,12 @@ real function get_entropy(rho,pres,mu_in,ieos) case default cgss = entropy(cgsrho,cgspres,mu_in,1) end select +! write(*,'(A,ES15.7)') 'ALI: 1 GET_ENTROPY: cgsS = ',cgss cgss = cgss/kboltz ! s/kb +! write(*,'(A,ES15.7)') 'ALI: 2 GET_ENTROPY: cgsS/kb = ',cgss + get_entropy = cgss/unit_ergg ! units in erg/grK, here it turns to code units +! write(*,'(A,ES15.7)') 'ALI: 3 GET_ENTROPY/kb/unit_ergg code units= ',get_entropy end function get_entropy From a60d31ae9aa1610347face09c2cc255221bd0e6d Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Mon, 6 Jul 2026 19:08:49 +1000 Subject: [PATCH 24/42] entropy is already to the power10 in get entropy in eos, changed that and now things work! I am testing it now! --- src/main/cons2primsolver.f90 | 27 --------------------------- src/main/eos.f90 | 16 ---------------- 2 files changed, 43 deletions(-) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 1719a5477..617e5f703 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -78,7 +78,6 @@ subroutine primitive2conservative(x,metrici,v,dens,u,P,rho,pmom,en,ien_type) use metric_tools, only:unpack_metric use io, only:error use eos, only:gmw,get_entropy - use units, only:unit_density,unit_ergg ! ali added real, intent(in) :: x(1:3),metrici(:,:,:) real, intent(in) :: dens,v(1:3),u,P real, intent(out) :: rho,pmom(1:3),en @@ -114,12 +113,7 @@ subroutine primitive2conservative(x,metrici,v,dens,u,P,rho,pmom,en,ien_type) if (ien_type == ien_etotal) then en = U0*enth*gvv + (1.+u)/U0 elseif (ien_type == ien_entropy_s) then - ! write(*,'(A,ES15.7)') 'before get_entropy in cons2primsolver:entropy (code units): ', en en = get_entropy(dens,P,gmw,ieos) - ! write(*,'(A,ES15.7)') 'after get_entropy in cons2primsolver:entropy (code units): ', en - ! write(*,'(A,ES15.7)') 'get_entropy:rho (cgs): ', dens*unit_density - ! write(*,'(A,ES15.7)') 'get_entropy:P (cgs): ', P*unit_ergg/unit_density - ! write(*,'(A,ES15.7)') 'get_entropy:u (cgs): ', u*unit_ergg else if (u > 0) then gam1 = 1. + P/(dens*u) @@ -163,12 +157,6 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie logical :: converged,have_eos_cache real :: gcov(0:3,0:3) ierr = 0 -! write(*,*) 'L158: ALI cons2primsolver: DEBUG ieos=', ieos, ' ien_type=', ien_type -! write(*,'(A,ES15.7)') 'L158:rho (cgs): ', dens*unit_density -! write(*,'(A,ES15.7)') 'L158:entropy (cgs): ', en*unit_ergg -! write(*,'(A,ES15.7)') 'L158:u (cgs): ', u*unit_ergg -! write(*,'(A,ES15.7)') 'L158:P (cgs): ', P*unit_ergg/unit_density -! write(*,'(A,ES15.7)') 'L158:Temperature (K): ', temp ! Get metric components from metric array call unpack_metric(metrici,gcov=gcov,gammaijUP=gammaijUP,alpha=alpha,betadown=betadown,betaUP=betaUP) @@ -205,19 +193,10 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie p = (gamma-1.)*dens*polyk elseif (ien_type == ien_entropy_s) then call get_p_from_rho_s(ieos,en,dens,gmw,P,temp) - ! write(*,*) 'ALI, en = ',en,' dens = ',dens,' P = ',P,' temp = ',temp - ! write(*,*) 'ENTER entropy branch' - ! write(*,*) 'ieos=', ieos select case(ieos) case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) - ! write(*,*) 'L200: ALI cons2primsolver, went into case 10 EOS MESA tables' - ! write(*,'(A,ES15.7)') 'L200:rho (cgs): ', dens*unit_density - ! write(*,'(A,ES15.7)') 'L200:entropy (cgs): ', en*unit_ergg - ! write(*,'(A,ES15.7)') 'L200:u (cgs): ', u*unit_ergg - ! write(*,'(A,ES15.7)') 'L200:P (cgs): ', P*unit_ergg/unit_density - ! write(*,'(A,ES15.7)') 'L200:Temperature (K): ', temp case (12) cgsdens = dens * unit_density @@ -285,12 +264,6 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie case (10) ! inputs and outputs are all in code units call get_u_from_rho_s(ieos,en,dens,u) - ! write(*,*) 'L263: ALI cons2primsolver, went into case 10 EOS MESA tables' - ! write(*,'(A,ES15.7)') 'L263:rho (cgs): ', dens*unit_density - ! write(*,'(A,ES15.7)') 'L263:entropy (cgs): ', en*unit_ergg - ! write(*,'(A,ES15.7)') 'L263:u (cgs): ', u*unit_ergg - ! write(*,'(A,ES15.7)') 'L263:P (cgs): ', P*unit_ergg/unit_density - ! write(*,'(A,ES15.7)') 'L263:Temperature (K): ', temp case (12) cgsdens = dens * unit_density diff --git a/src/main/eos.f90 b/src/main/eos.f90 index b3d79c0eb..096a64152 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -334,18 +334,6 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam if (present(gamma_local)) gamma_local = gam1 ! gamma is an output if (present(mu_local)) mu_local = 1./get_eos_1overmu_mesa(cgsrhoi,cgseni) if (ierr /= 0) call warning('eos_mesa','extrapolating off tables') - if (ierr /= 0) then - write(*,'(A)') '*** MESA EOS extrapolating off tables ***' - ! write(*,'(A,ES15.7)') 'rho (code units): ', rhoi - write(*,'(A,ES15.7)') 'rho (cgs): ', cgsrhoi - ! write(*,'(A,ES15.7)') 'u (code units): ', eni - write(*,'(A,ES15.7)') 'u (cgs): ', cgseni - ! write(*,'(A,ES15.7)') 'P (code units): ', presi - write(*,'(A,ES15.7)') 'P (cgs): ', cgspresi - write(*,'(A,F12.6)') 'gamma1: ', gam1 - write(*,'(A,ES15.7)') 'Temperature (K): ', temperaturei - call warning('eos_mesa','extrapolating off tables') - endif case(11) ! @@ -1137,7 +1125,6 @@ real function get_entropy(rho,pres,mu_in,ieos) real, intent(in) :: rho,pres,mu_in real :: cgsrho,cgspres,cgss -! write(*,*) 'ALI: GET_ENTROPY Called: ieos = ',ieos cgsrho = rho * unit_density cgspres = pres * unit_pressure select case (ieos) @@ -1148,12 +1135,9 @@ real function get_entropy(rho,pres,mu_in,ieos) case default cgss = entropy(cgsrho,cgspres,mu_in,1) end select -! write(*,'(A,ES15.7)') 'ALI: 1 GET_ENTROPY: cgsS = ',cgss cgss = cgss/kboltz ! s/kb -! write(*,'(A,ES15.7)') 'ALI: 2 GET_ENTROPY: cgsS/kb = ',cgss get_entropy = cgss/unit_ergg ! units in erg/grK, here it turns to code units -! write(*,'(A,ES15.7)') 'ALI: 3 GET_ENTROPY/kb/unit_ergg code units= ',get_entropy end function get_entropy From d39b37fe35f6dc406104da621d4945c852bba599 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 12:39:12 +1000 Subject: [PATCH 25/42] made reading gr tables only happen for the GR case --- src/main/eos_mesa.f90 | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index e0604f000..dddd06c9b 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -18,7 +18,7 @@ module eos_mesa ! use mesa_microphysics - + use dim, only:gr implicit none logical, private :: mesa_initialised = .false. @@ -57,14 +57,17 @@ subroutine init_eos_mesa(x,z,ierr) call get_eos_constants_mesa(ierr) if (ierr /= 0) return - mesa_eos_gr_prefix="output_rhos_" - call get_eos_constants_mesa_gr(ierr) + !!! 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) + call read_eos_mesa_gr(x,z,ierr) ! read GR tables + end if if (ierr /= 0) return call read_eos_mesa(x,z,ierr) call get_opacity_constants_mesa call read_opacity_mesa(x,z) - call read_eos_mesa_gr(x,z,ierr) ! read GR tables end subroutine init_eos_mesa From e04a696f9715b60b9f5f3234343b0237a46c1861 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 13:06:02 +1000 Subject: [PATCH 26/42] added a hidden variable called mesa_use_legacy_tables which is set to false by default to use thomas reichardt's old tables --- src/main/datafiles.f90 | 9 +++++++-- src/main/options.f90 | 5 ++++- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index 7cfdd0169..6d246ba25 100644 --- a/src/main/datafiles.f90 +++ b/src/main/datafiles.f90 @@ -56,13 +56,18 @@ end function find_phantom_datafile !+ !---------------------------------------------------------------- function map_dir_to_web(search_dir) result(url) + use options, only: mesa_use_legacy_tables character(len=*), intent(in) :: search_dir character(len=120) :: url !print*,' search_dir=',trim(search_dir) select case(search_dir) - case('data/eos/mesa') - url = 'https://zenodo.org/records/13148447/files/' + case('data/eos/mesa') !!! here this is if someone wants to cehck the older version of the mesa tables (Reichardt et al 2020), which are no longer used by default + if (mesa_use_legacy_tables) then + url = 'https://zenodo.org/records/13148447/files/' + else + url = 'https://zenodo.org/records/newlink' !! not uploaded yet + endif case('data/eos/shen') url = 'https://zenodo.org/records/13163155/files/' case('data/eos/helmholtz') diff --git a/src/main/options.f90 b/src/main/options.f90 index 4bc19e35d..94cdfdae7 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -49,7 +49,8 @@ module options ! library use logical, public :: write_files - +! MESA EOS table version (hidden variable, set via input file) + logical, public :: mesa_use_legacy_tables public :: set_default_options,write_options_output,read_options_output ! options from lower-level modules that can also be imported via options module @@ -125,6 +126,8 @@ subroutine set_default_options ! enable/disable writing output files write_files = .true. + ! use legacy MESA EOS tables + mesa_use_legacy_tables = .false. end subroutine set_default_options !----------------------------------------------------------------------- From 1ff204dce95723918cf693180961ec92650052d4 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 14:15:07 +1000 Subject: [PATCH 27/42] added a different link for EOS EMSA gr tables to be read, and also changed opacs table to be read separately for eos mesa as I have not updated them so they would just be downloaded from Thomas old link --- src/main/datafiles.f90 | 14 +++++++++----- src/main/eos_mesa_microphysics.f90 | 8 ++++---- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index 6d246ba25..692e183fe 100644 --- a/src/main/datafiles.f90 +++ b/src/main/datafiles.f90 @@ -63,11 +63,15 @@ function map_dir_to_web(search_dir) result(url) !print*,' search_dir=',trim(search_dir) select case(search_dir) case('data/eos/mesa') !!! here this is if someone wants to cehck the older version of the mesa tables (Reichardt et al 2020), which are no longer used by default - if (mesa_use_legacy_tables) then - url = 'https://zenodo.org/records/13148447/files/' - else - url = 'https://zenodo.org/records/newlink' !! not uploaded yet - endif + if (mesa_use_legacy_tables) then + url = 'https://zenodo.org/records/13148447/files/' + else + url = 'https://zenodo.org/records/newlink' !! not uploaded yet + endif + case('data/eos/mesa_entropy_gr') + url = 'https://zenodo.org/records/newlink' !! not uploaded yet + case('data/eos/mesa_opac') !!! same link as old eos/mesa, but this is only for downloading the opacity tables + url = 'https://zenodo.org/records/13148447/files/' case('data/eos/shen') url = 'https://zenodo.org/records/13163155/files/' case('data/eos/helmholtz') diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index ea2ba2c77..ab5b03a7a 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -80,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') @@ -114,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 @@ -364,7 +364,7 @@ subroutine get_eos_constants_mesa_gr(ierr) ! 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') + filename = find_phantom_datafile(filename,'eos/mesa_entropy_gr') ! Read constants from the header of first EoS tables open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted',iostat=ierr) @@ -521,7 +521,7 @@ subroutine read_eos_mesa_gr(x,z,ierr) !!!!!!!!!! have to rename these variables to not conflict with the above subroutine, but I can't be bothered to do that rn, so just ignore the fact that they have the same names as the above subroutine !!!!!!!!!! ! Find the EoS tables filename = trim(mesa_eos_gr_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' - filename = find_phantom_datafile(filename,'eos/mesa') + filename = find_phantom_datafile(filename,'eos/mesa_entropy_gr') ! Note that the data exists mesa_eos_gr_data_exists(i,j)=1 From 1da118edb0280d2ed40614dde119d664f1f9bbe7 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 14:21:03 +1000 Subject: [PATCH 28/42] minor addition --- src/main/datafiles.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index 692e183fe..ee9463b97 100644 --- a/src/main/datafiles.f90 +++ b/src/main/datafiles.f90 @@ -70,7 +70,7 @@ function map_dir_to_web(search_dir) result(url) endif case('data/eos/mesa_entropy_gr') url = 'https://zenodo.org/records/newlink' !! not uploaded yet - case('data/eos/mesa_opac') !!! same link as old eos/mesa, but this is only for downloading the opacity tables + 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/' From 2d019befd4c6d8e569e227e30dcddad877a77254 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 15:06:17 +1000 Subject: [PATCH 29/42] added some final touches. changed the logical to an integer for the eosmesa_version. also the gr tables will be uploaded in the same place that the classic tables will be --- src/main/datafiles.f90 | 20 ++++++++++++-------- src/main/eos.f90 | 4 ++-- src/main/eos_mesa.f90 | 2 +- src/main/eos_mesa_microphysics.f90 | 5 ++--- src/main/options.f90 | 6 +++--- 5 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index ee9463b97..be7676696 100644 --- a/src/main/datafiles.f90 +++ b/src/main/datafiles.f90 @@ -56,20 +56,24 @@ end function find_phantom_datafile !+ !---------------------------------------------------------------- function map_dir_to_web(search_dir) result(url) - use options, only: mesa_use_legacy_tables + use options, only: eosmesa_version character(len=*), intent(in) :: search_dir character(len=120) :: url !print*,' search_dir=',trim(search_dir) select case(search_dir) - case('data/eos/mesa') !!! here this is if someone wants to cehck the older version of the mesa tables (Reichardt et al 2020), which are no longer used by default - if (mesa_use_legacy_tables) then + 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/' - else - url = 'https://zenodo.org/records/newlink' !! not uploaded yet - endif - case('data/eos/mesa_entropy_gr') - url = 'https://zenodo.org/records/newlink' !! not uploaded yet + case(2) + url = 'https://zenodo.org/records/newlink/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') diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 096a64152..bff12845d 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1222,7 +1222,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.) cgspres = cgsrho*Rg*temp / mu case(10) - !!! For GR, not tested yet! + !!! For GR call get_eos_ptemp_from_rhos_mesa_gr(cgsrho,cgss,cgsP,temp) case (12) call get_idealplusrad_tempfromrhoS(cgsrho,cgss,mu,temp,cgspres,niter_out) @@ -1261,7 +1261,7 @@ subroutine get_u_from_rho_s(ieos,S,rho,u) select case (ieos) case(10) - !!! For GR, not tested yet! + !!! For GR call get_eos_u_from_rhos_mesa_gr(cgsrho,cgss,cgsu) case default diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index dddd06c9b..1cae542d7 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -14,7 +14,7 @@ module eos_mesa ! ! :Runtime parameters: None ! -! :Dependencies: mesa_microphysics, physcon +! :Dependencies: mesa_microphysics, physcon, dim ! use mesa_microphysics diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index ab5b03a7a..bceda10ea 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -364,7 +364,7 @@ subroutine get_eos_constants_mesa_gr(ierr) ! 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_entropy_gr') + 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) @@ -518,10 +518,9 @@ subroutine read_eos_mesa_gr(x,z,ierr) do j=1,mesa_eos_gr_nh write (hh,fmt1) mesa_eos_gr_h(j) ! X value as string -!!!!!!!!!! have to rename these variables to not conflict with the above subroutine, but I can't be bothered to do that rn, so just ignore the fact that they have the same names as the above subroutine !!!!!!!!!! ! Find the EoS tables filename = trim(mesa_eos_gr_prefix)//'z'//trim(zz)//'x'//trim(hh)//'.bindata' - filename = find_phantom_datafile(filename,'eos/mesa_entropy_gr') + filename = find_phantom_datafile(filename,'eos/mesa') ! Note that the data exists mesa_eos_gr_data_exists(i,j)=1 diff --git a/src/main/options.f90 b/src/main/options.f90 index 94cdfdae7..7813ea459 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -50,7 +50,7 @@ module options ! library use logical, public :: write_files ! MESA EOS table version (hidden variable, set via input file) - logical, public :: mesa_use_legacy_tables + integer, public :: eosmesa_version public :: set_default_options,write_options_output,read_options_output ! options from lower-level modules that can also be imported via options module @@ -126,8 +126,8 @@ subroutine set_default_options ! enable/disable writing output files write_files = .true. - ! use legacy MESA EOS tables - mesa_use_legacy_tables = .false. + ! 0 or 1 uses the older version of the mesa tables (Reichardt et al 2020), which are no longer used by default + eosmesa_version = 2 end subroutine set_default_options !----------------------------------------------------------------------- From a9709d124b7371ba6edcf56efd12e55834764de3 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 16:51:51 +1000 Subject: [PATCH 30/42] added zenodo link --- src/main/datafiles.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index be7676696..4da494e00 100644 --- a/src/main/datafiles.f90 +++ b/src/main/datafiles.f90 @@ -70,7 +70,7 @@ function map_dir_to_web(search_dir) result(url) case(0,1) url = 'https://zenodo.org/records/13148447/files/' case(2) - url = 'https://zenodo.org/records/newlink/files/' + url = 'https://zenodo.org/records/21253459/files/' case default stop 'Unknown eosmesa_version' end select From 7961b191fa1cec75611d4c3de5faac5270b77163 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 17:36:13 +1000 Subject: [PATCH 31/42] checkpoint --- src/main/datafiles.f90 | 3 ++- src/main/eos.f90 | 6 ++---- src/main/options.f90 | 5 +---- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/main/datafiles.f90 b/src/main/datafiles.f90 index 4da494e00..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 @@ -56,7 +58,6 @@ end function find_phantom_datafile !+ !---------------------------------------------------------------- function map_dir_to_web(search_dir) result(url) - use options, only: eosmesa_version character(len=*), intent(in) :: search_dir character(len=120) :: url diff --git a/src/main/eos.f90 b/src/main/eos.f90 index bff12845d..a42996d92 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -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 @@ -1208,7 +1208,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) real, intent(inout) :: temp real, intent(out) :: P integer, intent(in) :: ieos - real :: corr,df,f,temp_new,cgsrho,cgsp,cgss + real :: corr,df,f,temp_new,cgsrho,cgspres,cgss real, parameter :: eoserr=1e-12 integer :: niter integer, parameter :: nitermax = 1000 @@ -1248,7 +1248,6 @@ end subroutine get_p_from_rho_s subroutine get_u_from_rho_s(ieos,S,rho,u) use io, only:fatal use units, only:unit_density,unit_pressure,unit_ergg - use physcon, only:kboltz,avogadro use eos_mesa, only: get_eos_u_from_rhos_mesa_gr real, intent(in) :: S,rho real, intent(out) :: u @@ -1935,7 +1934,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/options.f90 b/src/main/options.f90 index 7813ea459..4bc19e35d 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -49,8 +49,7 @@ module options ! library use logical, public :: write_files -! MESA EOS table version (hidden variable, set via input file) - integer, public :: eosmesa_version + public :: set_default_options,write_options_output,read_options_output ! options from lower-level modules that can also be imported via options module @@ -126,8 +125,6 @@ subroutine set_default_options ! enable/disable writing output files write_files = .true. - ! 0 or 1 uses the older version of the mesa tables (Reichardt et al 2020), which are no longer used by default - eosmesa_version = 2 end subroutine set_default_options !----------------------------------------------------------------------- From 0b2a3a25ef6a191624c44d02297d03301ba2ca55 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 17:44:38 +1000 Subject: [PATCH 32/42] checkpoint --- src/main/eos.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index a42996d92..3d6afa099 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1212,10 +1212,12 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) real, parameter :: eoserr=1e-12 integer :: niter integer, parameter :: nitermax = 1000 + integer, intent(out), optional :: niter_out ! change to cgs unit cgsrho = rho*unit_density cgss = S*unit_ergg + if (present(niter_out)) niter_out = 0 select case (ieos) case (2,5) @@ -1223,7 +1225,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) cgspres = cgsrho*Rg*temp / mu case(10) !!! For GR - call get_eos_ptemp_from_rhos_mesa_gr(cgsrho,cgss,cgsP,temp) + 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 From 8d7e7042155f24506ee649407d4354b9d62f4438 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 19:04:40 +1000 Subject: [PATCH 33/42] missing declaration --- src/main/eos.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 3d6afa099..3ab29a854 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1201,7 +1201,8 @@ end subroutine get_rho_from_p_s 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_idealgasplusrad_tempfrompres,get_idealplusrad_pres + 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 From b288c9c7d68decbd21c43982b04618464d4584fb Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 19:26:19 +1000 Subject: [PATCH 34/42] something was deleted --- src/main/eos.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 3ab29a854..3453eb6d8 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -164,9 +164,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,10,12 or 11',& var='eos_type',val=real(eos_type)) endif From a55be3ebea45a859d3d50dd5241a8dd50598cdd5 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 20:09:39 +1000 Subject: [PATCH 35/42] removed useless variables --- src/main/eos.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 3453eb6d8..280112f28 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -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 @@ -1209,9 +1209,8 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp,niter_out) real, intent(inout) :: temp real, intent(out) :: P integer, intent(in) :: ieos - real :: corr,df,f,temp_new,cgsrho,cgspres,cgss + real :: cgsrho,cgspres,cgss real, parameter :: eoserr=1e-12 - integer :: niter integer, parameter :: nitermax = 1000 integer, intent(out), optional :: niter_out From 1365a7a08dac4c6fd9f2ee88e8c91adee91fbac4 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 22:10:52 +1000 Subject: [PATCH 36/42] minor debug --- src/main/eos.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 280112f28..747c0e86b 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1249,12 +1249,12 @@ end subroutine get_p_from_rho_s !----------------------------------------------------------------------- subroutine get_u_from_rho_s(ieos,S,rho,u) use io, only:fatal - use units, only:unit_density,unit_pressure,unit_ergg + 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,cgsp,cgss, cgsu + real :: cgsrho,cgss, cgsu ! change to cgs unit cgsrho = rho*unit_density From 4149beaf9091c094e9d13aad824ef616a10259d6 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 22:13:08 +1000 Subject: [PATCH 37/42] removed useless routines --- src/main/eos.f90 | 36 +----------------------------------- 1 file changed, 1 insertion(+), 35 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 747c0e86b..ca50924e8 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -69,7 +69,7 @@ module eos 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,get_u_from_rho_s - public :: calc_rho_from_PT,get_entropy,get_p_from_rho_s,get_entropy_vec, get_u_from_rho_s_vec + 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 public :: write_headeropts_eos,read_headeropts_eos @@ -1141,27 +1141,6 @@ real function get_entropy(rho,pres,mu_in,ieos) end function get_entropy -! vectorised version of get_entropy for use in relax_star -subroutine get_entropy_vec(rho, pres, mu_in, ieos, S) - use io, only: fatal - integer, intent(in) :: ieos - real, intent(in) :: rho(:), pres(:) - real, intent(in) :: mu_in - real, intent(out) :: S(:) - - integer :: i, n - - n = size(rho) - - if (size(pres) /= n .or. size(S) /= n) then - call fatal('get_entropy_vec','array size mismatch') - endif - - do i = 1, n - S(i) = get_entropy(rho(i), pres(i), mu_in, ieos) - end do - -end subroutine get_entropy_vec !----------------------------------------------------------------------- !+ ! Calculate density given pressure and entropy using Newton-Raphson @@ -1275,19 +1254,6 @@ subroutine get_u_from_rho_s(ieos,S,rho,u) end subroutine get_u_from_rho_s -! vectorised version of get_u_from_rho_s for use in relax_star -subroutine get_u_from_rho_s_vec(ieos,S,rho,u,n) - integer, intent(in) :: ieos, n - real, intent(in) :: S(n), rho(n) - real, intent(out) :: u(n) - - integer :: i - - do i = 1, n - call get_u_from_rho_s(ieos, S(i), rho(i), u(i)) - end do -end subroutine get_u_from_rho_s_vec - !----------------------------------------------------------------------- !+ ! Calculate mean molecular weight from X and Z, assuming complete From 511869fb25858a87b5f0b2d02a1ab760090d75cb Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 22:26:25 +1000 Subject: [PATCH 38/42] changed some issues --- src/main/eos.f90 | 2 +- src/main/eos_mesa.f90 | 14 +++++++++----- src/main/eos_mesa_microphysics.f90 | 6 ++++-- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index ca50924e8..923ab3653 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -166,7 +166,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam ! 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,10,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 diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index 1cae542d7..aa0773a4c 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -61,11 +61,15 @@ subroutine init_eos_mesa(x,z,ierr) !!! 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) - call read_eos_mesa_gr(x,z,ierr) ! read GR tables + 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 - if (ierr /= 0) return + call read_eos_mesa(x,z,ierr) + if (ierr /= 0) return call get_opacity_constants_mesa call read_opacity_mesa(x,z) @@ -141,7 +145,7 @@ subroutine get_eos_pressure_temp_mesa(den,eint,pres,temp) end subroutine get_eos_pressure_temp_mesa -!---------------------------------------------------------------- NOT TESTED YET! +!---------------------------------------------------------------- !+ ! subroutine returns pressure and temperature as ! a function of density/entropy for GR tables @@ -156,7 +160,7 @@ subroutine get_eos_ptemp_from_rhos_mesa_gr(den,s,pres,temp) end subroutine get_eos_ptemp_from_rhos_mesa_gr -!---------------------------------------------------------------- NOT TESTED YET! +!---------------------------------------------------------------- !+ ! subroutine returns internal energy as ! a function of density/entropy for GR tables diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index bceda10ea..12093a705 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -368,7 +368,7 @@ subroutine get_eos_constants_mesa_gr(ierr) ! 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 (again, I wonder if it's best to make another subroutine for this) +! allocate GR tables if (ierr /= 0) return read(fnum) mesa_eos_gr_ns, mesa_eos_gr_nrho, mesa_eos_gr_nvar2 close(fnum) @@ -847,7 +847,9 @@ 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 (I wonder if it's best to make another subroutine for deallocating the GR tables) +!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) From b50ffd7cc1235ad0455158ecfc83aad3b7b3d9af Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Wed, 8 Jul 2026 23:14:36 +1000 Subject: [PATCH 39/42] fixed important bug in eos_mesa on extrapolation --- src/main/eos.f90 | 4 ++-- src/main/eos_mesa_microphysics.f90 | 25 ++++++++++++++----------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 923ab3653..64cb8c54b 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1101,8 +1101,8 @@ 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 (output is not logs, - ! because for vars >5 in the tables the log becomes the actual value) + ! 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,entropy,ierr) else diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index 12093a705..e2bdd0caf 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -629,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 ! this is applied only for the first four columns of the data which are stored in log10 form. + 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 @@ -679,16 +681,17 @@ pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) 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 (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 (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 = 10.d0**((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)) + 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 From fc05a5ecf45b9c64730c74566840f8062e1f3c24 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Thu, 16 Jul 2026 19:27:21 +1000 Subject: [PATCH 40/42] minor bug and merge --- src/main/cons2primsolver.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 617e5f703..65d7c1569 100644 --- a/src/main/cons2primsolver.f90 +++ b/src/main/cons2primsolver.f90 @@ -193,6 +193,10 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie p = (gamma-1.)*dens*polyk elseif (ien_type == ien_entropy_s) then call get_p_from_rho_s(ieos,en,dens,gmw,P,temp) + dens_eos = dens + p_eos = P + temp_eos = temp + have_eos_cache = .true. select case(ieos) case (10) ! inputs and outputs are all in code units From cc08c59062d88ffb43d61df3037c0ba8ee298789 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Fri, 17 Jul 2026 09:50:24 +1000 Subject: [PATCH 41/42] minor debugs, plus rerun prompt to see the ssl nasa jpl problem --- src/main/eos.f90 | 4 ++-- src/main/eos_mesa.f90 | 2 +- src/main/eos_mesa_microphysics.f90 | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index f03f9fa91..75d090d6c 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1117,14 +1117,14 @@ 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 +! 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 integer, intent(in) :: ieos real, intent(in) :: rho,pres,mu_in real :: cgsrho,cgspres,cgss - + cgsrho = rho * unit_density cgspres = pres * unit_pressure select case (ieos) diff --git a/src/main/eos_mesa.f90 b/src/main/eos_mesa.f90 index aa0773a4c..066ac80eb 100644 --- a/src/main/eos_mesa.f90 +++ b/src/main/eos_mesa.f90 @@ -160,7 +160,7 @@ subroutine get_eos_ptemp_from_rhos_mesa_gr(den,s,pres,temp) end subroutine get_eos_ptemp_from_rhos_mesa_gr -!---------------------------------------------------------------- +!---------------------------------------------------------------- !+ ! subroutine returns internal energy as ! a function of density/entropy for GR tables diff --git a/src/main/eos_mesa_microphysics.f90 b/src/main/eos_mesa_microphysics.f90 index e2bdd0caf..a2931c9bc 100644 --- a/src/main/eos_mesa_microphysics.f90 +++ b/src/main/eos_mesa_microphysics.f90 @@ -52,7 +52,7 @@ module mesa_microphysics integer, allocatable :: mesa_eos_gr_data_exists(:,:) real, allocatable :: mesa_gr_ds_data(:,:,:) real, allocatable :: mesa_eos_gr0(:,:,:,:) - real, allocatable :: mesa_gr_ds_data0(:,:,:,:,:) + real, allocatable :: mesa_gr_ds_data0(:,:,:,:,:) public :: get_opacity_constants_mesa public :: read_opacity_mesa @@ -368,7 +368,7 @@ subroutine get_eos_constants_mesa_gr(ierr) ! 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 +! allocate GR tables if (ierr /= 0) return read(fnum) mesa_eos_gr_ns, mesa_eos_gr_nrho, mesa_eos_gr_nvar2 close(fnum) @@ -681,7 +681,7 @@ pure subroutine getvalue_mesa_gr(rho,s,ivout,vout,ierr) 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 + ! 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. From 840dd11df8e64e4b8c4713169c19e3b6fe292750 Mon Sep 17 00:00:00 2001 From: Ali Pourmand Date: Fri, 17 Jul 2026 17:18:47 +1000 Subject: [PATCH 42/42] added readme in mesa_opac --- data/eos/mesa_opac/README | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 data/eos/mesa_opac/README 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