#include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" ! eoskey: ! 1 --- polytropic EOS ! 2 --- gamma-law EOS ! 3 --- hybrid EOS ! 4 --- finite-T microphysical NSE EOS subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,& rho,eps,temp,ye,press,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,keytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) CCTK_REAL, intent(out) :: press(npoints) ! local vars integer :: i character(256) :: warnstring real*8 :: hybrid_local_gamma, hybrid_local_k, & hybrid_p_poly, hybrid_p_th real*8,parameter :: zero = 0.0d0 ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt,xdpderho,xdpdrhoe ! temporary vars for coldeos + gamma law integer :: ir real*8 :: press_cold, press_th real*8 :: eps_cold, eps_th real*8 :: gamma anyerr = 0 keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = poly_k * & rho(i)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints press(i) = poly_k * & rho(i)**poly_gamma enddo case (2) ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = gl_k * & rho(i)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i) enddo case (3) ! hybrid EOS do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 hybrid_local_k = hybrid_k2 else hybrid_local_gamma = hybrid_gamma1 hybrid_local_k = hybrid_k1 endif hybrid_p_poly = hybrid_local_k * & rho(i)**hybrid_local_gamma hybrid_p_th = - hybrid_local_k * (hybrid_gamma_th - 1.d0) / & (hybrid_local_gamma - 1.0d0) * rho(i)**hybrid_local_gamma + & (hybrid_gamma_th - 1.0d0) * rho(i) * eps(i) - & (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & hybrid_k1 * & hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) * rho(i) hybrid_p_th = max(zero, hybrid_p_th) press(i) = hybrid_p_poly + hybrid_p_th enddo case (4) ! nuc eos do i=1,npoints xrho = rho(i) * inv_rho_gf xtemp = temp(i) xye = ye(i) xenr = eps(i) * inv_eps_gf call nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& keytemp,keyerr(i),rf_precision) if(keyerr(i).ne.0) then anyerr = 1 endif if(keytemp.eq.1) then eps(i) = xenr * eps_gf else temp(i) = xtemp endif press(i) = xprs * press_gf enddo case (5) ! cold tabular EOS with gamma law do i=1,npoints if(rho(i).lt.coldeos_rhomin) then press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i) cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 else xrho = log10(rho(i)) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1) eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1) if(keytemp.eq.0) then eps_th = max(0.0d0,eps(i) - eps_cold) else if (keytemp.eq.1) then eps_th = 0.0d0 eps(i) = eps_cold endif press_cold = coldeos_kappa * rho(i)**gamma press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th press(i) = press_cold + press_th enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_Press subroutine EOS_Omni_EOS_PressOMP(eoskey,keytemp,rf_precision,npoints,& rho,eps,temp,ye,press,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,keytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) CCTK_REAL, intent(out) :: press(npoints) ! local vars integer :: i character(256) :: warnstring real*8 :: hybrid_local_gamma, hybrid_local_k, & hybrid_p_poly, hybrid_p_th real*8,parameter :: zero = 0.0d0 ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt,xdpderho,xdpdrhoe ! temporary vars for coldeos + gamma law integer :: ir real*8 :: press_cold, press_th real*8 :: eps_cold, eps_th real*8 :: gamma anyerr = 0 keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS if(keytemp.eq.1) then !$OMP PARALLEL DO do i=1,npoints eps(i) = poly_k * & rho(i)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo !$OMP END PARALLEL DO endif !$OMP PARALLEL DO do i=1,npoints press(i) = poly_k * & rho(i)**poly_gamma enddo !$OMP END PARALLEL DO case (2) ! gamma-law EOS if(keytemp.eq.1) then !$OMP PARALLEL DO do i=1,npoints eps(i) = gl_k * & rho(i)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo !$OMP END PARALLEL DO endif !$OMP PARALLEL DO do i=1,npoints press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i) enddo !$OMP END PARALLEL DO case (3) ! hybrid EOS !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,& !$OMP hybrid_p_poly, hybrid_p_th) do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 hybrid_local_k = hybrid_k2 else hybrid_local_gamma = hybrid_gamma1 hybrid_local_k = hybrid_k1 endif hybrid_p_poly = hybrid_local_k * & rho(i)**hybrid_local_gamma hybrid_p_th = - hybrid_local_k * (hybrid_gamma_th - 1.d0) / & (hybrid_local_gamma - 1.0d0) * rho(i)**hybrid_local_gamma + & (hybrid_gamma_th - 1.0d0) * rho(i) * eps(i) - & (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & hybrid_k1 * & hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) * rho(i) hybrid_p_th = max(zero, hybrid_p_th) press(i) = hybrid_p_poly + hybrid_p_th enddo !$OMP END PARALLEL DO case (4) ! nuc eos !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs) do i=1,npoints xrho = rho(i) * inv_rho_gf xtemp = temp(i) xye = ye(i) xenr = eps(i) * inv_eps_gf call nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& keytemp,keyerr(i),rf_precision) if(keyerr(i).ne.0) then !$OMP CRITICAL anyerr = 1 !$OMP END CRITICAL endif if(keytemp.eq.1) then eps(i) = xenr * eps_gf else temp(i) = xtemp endif press(i) = xprs * press_gf enddo !$OMP END PARALLEL DO case (5) ! cold tabular EOS with gamma law !$OMP PARALLEL DO PRIVATE(xrho,ir,gamma,eps_cold,eps_th, & !$OMP press_cold,press_th) do i=1,npoints if(rho(i).lt.coldeos_rhomin) then press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i) cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 !$OMP CRITICAL anyerr = 1 !$OMP END CRITICAL else xrho = log10(rho(i)) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1) eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1) if(keytemp.eq.0) then eps_th = max(0.0d0,eps(i) - eps_cold) else if (keytemp.eq.1) then eps_th = 0.0d0 eps(i) = eps_cold endif press_cold = coldeos_kappa * rho(i)**gamma press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th press(i) = press_cold + press_th enddo !$OMP END PARALLEL DO case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_PressOMP subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,& rho,eps,temp,ye,dpdrhoe,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,keytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) CCTK_REAL, intent(out) :: dpdrhoe(npoints) ! local vars integer :: i character(256) :: warnstring real*8 :: hybrid_local_gamma, hybrid_local_k, & hybrid_dp_poly, hybrid_dp_th1, hybrid_dp_th2 real*8,parameter :: zero = 0.0d0 ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt,xdpderho,xdpdrhoe ! temporary vars for cold tabulated EOS + gamma law integer :: ir real*8 :: press_cold, press_th real*8 :: eps_cold, eps_th real*8 :: gamma, cs2, h anyerr = 0 keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = poly_k * & rho(i)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints dpdrhoe(i) = poly_k * & poly_gamma * & rho(i) ** (poly_gamma - 1.d0) enddo case (2) ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = gl_k * & rho(i)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints dpdrhoe(i) = (gl_gamma-1.0d0) * & eps(i) enddo case (3) ! hybrid EOS do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 hybrid_local_k = hybrid_k2 else hybrid_local_gamma = hybrid_gamma1 hybrid_local_k = hybrid_k1 endif hybrid_dp_poly = hybrid_local_gamma * & hybrid_local_k * rho(i)**(hybrid_local_gamma - 1.0d0) hybrid_dp_th1 = - hybrid_local_gamma * hybrid_local_k * & (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) * & rho(i)**(hybrid_local_gamma - 1.d0) hybrid_dp_th2 = (hybrid_gamma_th - 1.d0) * eps(i) & - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & hybrid_k1 * & hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) dpdrhoe(i) = hybrid_dp_poly + max(0.0d0,hybrid_dp_th1 + hybrid_dp_th2) enddo case (4) do i=1,npoints xrho = rho(i) * inv_rho_gf xtemp = temp(i) xye = ye(i) xenr = eps(i) * inv_eps_gf call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,& keytemp,keyerr(i),rf_precision) if(keyerr(i).ne.0) then anyerr = 1 endif if(keytemp.eq.1) then eps(i) = xenr * eps_gf else temp(i) = xtemp endif dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf enddo case (5) ! with the cold eos we have to assume P = P(rho), so ! by definition dPdrho is at constant internal energy ! and entropy (the latter, because T=0) do i=1,npoints if(rho(i).lt.coldeos_rhomin) then dpdrhoe(i) = coldeos_low_kappa * coldeos_low_gamma * & rho(i)**(coldeos_low_gamma-1.0d0) cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 else xrho = log10(rho(i)) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1) ! this is the cold speed of sound squared cs2 = (coldeos_cs2(ir) - coldeos_cs2(ir-1)) / & (coldeos_cs2(ir) - coldeos_cs2(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_cs2(ir-1) eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1) if(keytemp.eq.0) then eps_th = max(0.0d0,eps(i) - eps_cold) else if (keytemp.eq.1) then eps_th = 0.0d0 eps(i) = eps_cold endif press_cold = coldeos_kappa * rho(i)**gamma press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th ! this is the cold enthalpy, because it belongs to the ! cold speed of sound squared h = 1.0d0 + eps_cold + press_cold / rho(i) dpdrhoe(i) = cs2*h + press_th / rho(i) enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_DPressByDRho subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,& rho,eps,temp,ye,dpdepsrho,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,keytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) CCTK_REAL, intent(out) :: dpdepsrho(npoints) ! local vars integer :: i character(256) :: warnstring ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt,xdpderho,xdpdrhoe ! temporary vars for coldeos + gamma law integer :: ir real*8 :: eps_cold, eps_th anyerr = 0 keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = poly_k * & rho(i)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints dpdepsrho(i) = 0.0d0 enddo case (2) ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = gl_k * & rho(i)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints dpdepsrho(i) = (gl_gamma - 1.0d0) * & rho(i) enddo case (3) ! hybrid EOS do i=1,npoints dpdepsrho(i) = (hybrid_gamma_th - 1.0d0) * rho(i) enddo case (4) ! nuc_eos do i=1,npoints xrho = rho(i) * inv_rho_gf xtemp = temp(i) xye = ye(i) xenr = eps(i) * inv_eps_gf call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,& keytemp,keyerr(i),rf_precision) if(keyerr(i).ne.0) then anyerr = 1 endif if(keytemp.eq.1) then eps(i) = xenr * eps_gf else temp(i) = xtemp endif dpdepsrho(i) = xdpderho * press_gf * inv_eps_gf enddo case (5) ! with the cold eos we have to assume P = P(rho), so ! only the gamma law has non-zero dPdeps do i=1,npoints if(rho(i).lt.coldeos_rhomin) then dpdepsrho(i) = 0.0d0 cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 else xrho = log10(rho(i)) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1) if(keytemp.eq.0) then eps_th = max(0.0d0,eps(i) - eps_cold) else if (keytemp.eq.1) then eps_th = 0.0d0 eps(i) = eps_cold endif dpdepsrho(i) = coldeos_thfac * (coldeos_gammath - 1.0d0)*rho(i) enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_DPressByDEps subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,& rho,eps,temp,ye,cs2,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,keytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) CCTK_REAL, intent(out) :: cs2(npoints) ! local vars integer :: i character(256) :: warnstring real*8 :: xpress,xdpdrhoe,xdpderho real*8 :: hybrid_local_gamma, hybrid_local_k, & hybrid_p_poly, hybrid_p_th real*8,parameter :: zero = 0.0d0 ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt ! temporary vars for cold tabulated EOS + gamma law integer :: ir real*8 :: press_cold, press_th real*8 :: eps_cold, eps_th real*8 :: gamma, cs2_cold, cs2_th, h, h_cold anyerr = 0 keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = poly_k * & rho(i)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints xpress = poly_k * & rho(i)**(poly_gamma) cs2(i) = poly_gamma * xpress / rho(i) / & (1 + eps(i) + xpress/rho(i)) enddo case (2) ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints eps(i) = gl_k * & rho(i)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif do i=1,npoints xpress = (gl_gamma-1.0d0)*rho(i)*eps(i) xdpdrhoe = (gl_gamma-1.0d0)*eps(i) xdpderho = (gl_gamma-1.0d0)*rho(i) cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / & (1.0d0 + eps(i) + xpress/rho(i)) enddo case(3) ! hybrid EOS do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 hybrid_local_k = hybrid_k2 else hybrid_local_gamma = hybrid_gamma1 hybrid_local_k = hybrid_k1 endif ! first calculate the pressure hybrid_p_poly = hybrid_local_k * & rho(i)**hybrid_local_gamma hybrid_p_th = - hybrid_local_k * (hybrid_gamma_th - 1.d0) / & (hybrid_local_gamma - 1.0d0) * rho(i)**hybrid_local_gamma + & (hybrid_gamma_th - 1.0d0) * rho(i) * eps(i) - & (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & hybrid_k1 * & hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) * rho(i) hybrid_p_th = max(zero, hybrid_p_th) xpress = hybrid_p_poly + hybrid_p_th cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / & rho(i) / (1.0d0 + eps(i) + xpress/rho(i)) enddo case(4) ! nuc_eos do i=1,npoints xrho = rho(i) * inv_rho_gf xtemp = temp(i) xye = ye(i) xenr = eps(i) * inv_eps_gf call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,& keytemp,keyerr(i),rf_precision) if(keyerr(i).ne.0) then anyerr = 1 endif if(keytemp.eq.1) then eps(i) = xenr * eps_gf else temp(i) = xtemp endif cs2(i) = xcs2 * cliteinv2 / & (1.0d0 + eps(i) + xprs * press_gf / rho(i)) enddo case (5) ! with the cold eos we have to assume P = P(rho), so ! by definition dPdrho is at constant internal energy ! and entropy (the latter, because T=0) do i=1,npoints if(rho(i).lt.coldeos_rhomin) then xprs = coldeos_low_kappa * rho(i)**coldeos_low_gamma eps(i) = xprs / (coldeos_low_gamma - 1.0d0) / rho(i) cs2(i) = coldeos_low_kappa * coldeos_low_gamma * & rho(i)**(coldeos_low_gamma-1.0d0) / & (1.0d0 + eps(i) + xprs) cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 else xrho = log10(rho(i)) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1) ! this is the cold speed of sound squared cs2_cold = (coldeos_cs2(ir) - coldeos_cs2(ir-1)) / & (coldeos_cs2(ir) - coldeos_cs2(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_cs2(ir-1) eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1) if(keytemp.eq.0) then eps_th = max(0.0d0,eps(i) - eps_cold) else if (keytemp.eq.1) then eps_th = 0.0d0 eps(i) = eps_cold endif press_cold = coldeos_kappa * rho(i)**gamma press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th xdpdrhoe = coldeos_thfac*(coldeos_gammath - 1.0d0)*eps_th xdpderho = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i) cs2_th = (xdpdrhoe + press_th * xdpderho / (rho(i)**2)) h = 1.0d0 + eps(i) + (press_cold+press_th) / rho(i) h_cold = 1.0d0 + eps_cold + press_cold / rho(i) cs2(i) = (cs2_cold * h_cold + cs2_th) / h enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_cs2 subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,rf_precision,npoints,& rho,eps,temp,ye,press,xeps,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,keytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: rho(npoints),ye(npoints),press(npoints) CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) CCTK_REAL, intent(out) :: xeps(npoints) ! local vars integer :: i character(256) :: warnstring ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt,xdpderho,xdpdrhoe if(keytemp.eq.1) then anyerr = 1 keyerr(:) = -1 else anyerr = 0 keyerr(:) = 0 endif select case (eoskey) case (1) ! polytropic EOS do i=1,npoints xeps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i) enddo case (2) ! gamma-law EOS do i=1,npoints xeps(i) = press(i) / (gl_gamma - 1.0d0) / rho(i) enddo case (3) ! hybrid EOS write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for hybrid EOS" call CCTK_WARN(0,warnstring) case (4) ! nuc EOS write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for nuc_eos yet" call CCTK_WARN(0,warnstring) case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_eps_from_press subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,& ikeytemp,rf_precision,& npoints,rho,eps,temp,entropy,ye,press,keyerr,anyerr) ! This routine returns the density and spec. internal energy based on inputs of ! pressure, temperature, and Y_e. ! The current implementation is robust, but very slow; it should be used only ! for initial data where speed does not matter. use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS CCTK_INT, intent(in) :: eoskey,ikeytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision CCTK_REAL, intent(in) :: ye(npoints),press(npoints) CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints),eps(npoints) CCTK_REAL, intent(inout) :: entropy(npoints) ! local vars integer :: i character(256) :: warnstring integer :: keytemp ! temporary vars for nuc_eos real*8 :: xrho,xye,xtemp,xenr,xent real*8 :: xprs,xmunu,xcs2 real*8 :: xdedt,xdpderho,xdpdrhoe ! temporary vars for cold tabulated EOS integer :: ir real*8 :: gamma ! helper vars logical :: cont integer :: counter real*8 :: rho_guess, rho_guess2 real*8 :: press_guess,press_guess2, mydpdrho real*8 :: fac keytemp = ikeytemp anyerr = 0 keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS do i=1,npoints rho(i) = (press(i) / poly_k)**(1.0/poly_gamma) eps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i) enddo case (2) ! gamma-law EOS do i=1,npoints rho(i) = press(i) / ((gl_gamma - 1.0d0)*eps(i)) enddo case (3) ! hybrid EOS write(warnstring,*) "EOS_Omni_RestMassDensityFromPressEpsTemp not supported for hybrid EOS" call CCTK_WARN(0,warnstring) case (4) ! nuc EOS if(ikeytemp.eq.2) then call CCTK_WARN(0,"This function does not work yet when coming in with entropy") else if(ikeytemp.eq.1) then keytemp = 1 else call CCTK_WARN(0,"This function does not work yet when coming in with this keytemp") endif keytemp = 1 do i=1,npoints fac = 1.0d0 counter = 0 cont = .true. xprs = press(i) * inv_press_gf press_guess = xprs rho_guess = rho(i) * inv_rho_gf xprs = press(i) * inv_press_gf xtemp = temp(i) xye = ye(i) do while(cont) counter = counter + 1 rho_guess2 = rho_guess * 1.0001d0 call nuc_eos_short(rho_guess2,xtemp,xye,xenr,press_guess2,& xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,& keytemp,keyerr(i),rf_precision) call nuc_eos_short(rho_guess,xtemp,xye,xenr,press_guess,& xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,& keytemp,keyerr(i),rf_precision) mydpdrho = (press_guess2-press_guess)/(rho_guess2-rho_guess) if (mydpdrho.lt.0.0d0) then !$OMP CRITICAL write(warnstring,"(A25,1P10E15.6)") "Issue with table, dpdrho.lt.0",& rho_guess,xtemp,xye call CCTK_WARN(0,warnstring) !$OMP END CRITICAL endif if (abs(1.0d0-press_guess/xprs).lt.rf_precision) then cont = .false. rho(i) = rho_guess*rho_gf eps(i) = xenr*eps_gf else if (fac*(xprs-press_guess)/mydpdrho/rho_guess.gt.0.1d0) then rho_guess = 0.99d0*rho_guess else rho_guess = rho_guess + fac*(xprs-press_guess)/mydpdrho endif endif if (counter.gt.100) then fac = 0.01d0 endif if (rho_guess.lt.1.0d3.or.counter.gt.100000) then keyerr(i) = 473 anyerr = 1 return endif enddo enddo case (5) ! cold tabulate EOS ! deal only with case in which thermal pressure is zero if(keytemp.ne.1) then call CCTK_WARN(0,"finding rho(press) for tabulated cold EOS only possible if keytemp=1") endif do i=1,npoints fac = 1.0d0 counter = 0 cont = .true. xprs = press(i) press_guess = xprs rho_guess = rho(i) xprs = press(i) do while(cont) counter = counter + 1 rho_guess2 = rho_guess * 1.0001d0 ! press 2 if(rho_guess2.lt.coldeos_rhomin) then keyerr(i) = 104 anyerr = 1 else if(rho_guess2.gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 else xrho = log10(rho_guess2) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1) press_guess2 = coldeos_kappa * rho_guess2**gamma ! end press 2 ! press 1 if(rho_guess.lt.coldeos_rhomin) then keyerr(i) = 104 anyerr = 1 else if(rho_guess.gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 else xrho = log10(rho_guess) ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi ) endif ir = max(2, min(ir,coldeos_nrho)) gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1) press_guess = coldeos_kappa * rho_guess**gamma ! end press 1 ! derivative mydpdrho = (press_guess2-press_guess)/(rho_guess2-rho_guess) if (mydpdrho.lt.0.0d0) then !$OMP CRITICAL write(warnstring,"(A25,1P10E15.6)") "Issue with table, dpdrho.lt.0",& rho_guess,rho_guess2 call CCTK_WARN(0,warnstring) !$OMP END CRITICAL endif if (abs(1.0d0-press_guess/xprs).lt.rf_precision) then cont = .false. rho(i) = rho_guess eps(i) = (coldeos_eps(ir) - coldeos_eps(ir-1)) / & (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * & (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1) else if (fac*(xprs-press_guess)/mydpdrho/rho_guess.gt.0.1d0) then rho_guess = 0.99d0*rho_guess else rho_guess = rho_guess + fac*(xprs-press_guess)/mydpdrho endif endif if (counter.gt.100) then fac = 0.01d0 endif if (rho_guess.lt.1.0d-15.or.counter.gt.100000) then keyerr(i) = 473 anyerr = 1 return endif enddo enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select end subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt