diff options
Diffstat (limited to 'src/EOS_Omni_SingleVarCalls.F90')
| -rw-r--r-- | src/EOS_Omni_SingleVarCalls.F90 | 119 |
1 files changed, 106 insertions, 13 deletions
diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90 index aa0f79d..60c656d 100644 --- a/src/EOS_Omni_SingleVarCalls.F90 +++ b/src/EOS_Omni_SingleVarCalls.F90 @@ -24,8 +24,11 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,npoints,& CCTK_REAL, intent(out) :: press(npoints) ! local vars - integer :: i - character(256) :: warnstring + integer :: i + character(256) :: warnstring + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & + hybrid_p_poly, hybrid_p_th + real*8,parameter :: zero = 0.0d0 anyerr = 0 keyerr(:) = 0 @@ -57,6 +60,29 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,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_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + hybrid_p_poly = press_gf * hybrid_local_k_cgs * & + (rho(i) * inv_rho_gf)**hybrid_local_gamma + hybrid_p_th = - press_gf * hybrid_local_k_cgs * (hybrid_gamma_th - 1.d0) / & + (hybrid_local_gamma - 1.0d0) * (rho(i) * inv_rho_gf)**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) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + 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 DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -81,8 +107,11 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,npoints,& CCTK_REAL, intent(out) :: dpdrhoe(npoints) ! local vars - integer :: i - character(256) :: warnstring + integer :: i + character(256) :: warnstring + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & + hybrid_dp_poly, hybrid_dp_th1, hybrid_dp_th2 + real*8,parameter :: zero = 0.0d0 anyerr = 0 keyerr(:) = 0 @@ -115,6 +144,33 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,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_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + hybrid_dp_poly = hybrid_local_gamma * press_gf * & + hybrid_local_k_cgs * rho(i)**(hybrid_local_gamma - 1.0d0) * & + inv_rho_gf**hybrid_local_gamma + + hybrid_dp_th1 = - hybrid_local_gamma * press_gf * hybrid_local_k_cgs * & + (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) * & + rho(i)**(hybrid_local_gamma - 1.d0) * inv_rho_gf**hybrid_local_gamma + + 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) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) + + dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2 + + enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" @@ -138,9 +194,8 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,npoints,& CCTK_REAL, intent(out) :: dpdepsrho(npoints) ! local vars - integer :: i - character(256) :: warnstring - + integer :: i + character(256) :: warnstring anyerr = 0 keyerr(:) = 0 @@ -170,6 +225,11 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,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 DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" @@ -193,9 +253,12 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,npoints,& CCTK_REAL, intent(out) :: cs2(npoints) ! local vars - integer :: i - character(256) :: warnstring - real*8 :: xpress,xdpdrhoe,xdpderho + integer :: i + character(256) :: warnstring + real*8 :: xpress,xdpdrhoe,xdpderho + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & + hybrid_p_poly, hybrid_p_th + real*8,parameter :: zero = 0.0d0 anyerr = 0 keyerr(:) = 0 @@ -232,7 +295,31 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,npoints,& 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_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + ! first calculate the pressure + hybrid_p_poly = press_gf * hybrid_local_k_cgs * & + (rho(i) * inv_rho_gf)**hybrid_local_gamma + hybrid_p_th = - press_gf * hybrid_local_k_cgs * (hybrid_gamma_th - 1.d0) / & + (hybrid_local_gamma - 1.0d0) * (rho(i) * inv_rho_gf)**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) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + 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 DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -272,14 +359,16 @@ subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,npoints,& ! polytropic EOS do i=1,npoints xeps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i) - eps(i) = xeps(i) enddo case (2) ! gamma-law EOS do i=1,npoints xeps(i) = press(i) / (gl_gamma - 1.0d0) / rho(i) - eps(i) = xeps(i) enddo + case (3) + ! hybrid EOS + write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for hybrid EOS" + call CCTK_WARN(0,warnstring) case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -325,6 +414,10 @@ subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,npoints,& do i=1,npoints rho(i) = press(i) / ((gl_gamma - 1.0d0)*eps(i)) enddo + case (3) + ! hybrid EOS + write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress not supported for hybrid EOS" + call CCTK_WARN(0,warnstring) case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) |