[Commits] [svn:einsteintoolkit] incoming/EOS_Omni/ (Rev. 3)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Fri Jul 2 16:41:51 CDT 2010
User: cott
Date: 2010/07/02 04:41 PM
Modified:
/EOS_Omni/
param.ccl
/EOS_Omni/src/
EOS_Omni_Module.F90, EOS_Omni_SingleVarCalls.F90, EOS_Omni_Startup.F90
Log:
* add hybrid EOS
File Changes:
Directory: /EOS_Omni/
=====================
File [modified]: param.ccl
Delta lines: +33 -2
===================================================================
--- EOS_Omni/param.ccl 2010-07-02 21:41:36 UTC (rev 2)
+++ EOS_Omni/param.ccl 2010-07-02 21:41:51 UTC (rev 3)
@@ -4,18 +4,23 @@
restricted:
# poly EOS
-CCTK_REAL poly_gamma "Adiabatic Index for poly EOS"
+REAL poly_gamma "Adiabatic Index for poly EOS"
{
: :: ""
} 2.0
+REAL poly_gamma_ini "Initial Adiabatic Index for poly EOS"
+{
+ : :: ""
+} 2.0
+
REAL poly_k "Polytropic constant in c=G=Msun=1"
{
: :: ""
} 100.0
# gamma-law EOS
-CCTK_REAL gl_gamma "Adiabatic Index for gamma-law EOS"
+REAL gl_gamma "Adiabatic Index for gamma-law EOS"
{
: :: ""
} 2.0
@@ -24,3 +29,29 @@
{
: :: ""
} 100.0
+
+# hybrid EOS
+REAL hybrid_gamma1 "subnuclear adiabatic Index for hybrid EOS"
+{
+ : :: ""
+} 1.325
+
+REAL hybrid_gamma2 "subnuclear adiabatic Index for hybrid EOS"
+{
+ : :: ""
+} 2.5
+
+REAL hybrid_gamma_th "Thermal gamma for hybrid EOS"
+{
+ : :: ""
+} 1.5
+
+REAL hybrid_k1 "Polytropic constant in c=G=Msun=1 for hybrid EOS"
+{
+ : :: ""
+} 0.4640517
+
+REAL hybrid_rho_nuc "Density at which to switch betwen gammas; c=G=Msun=1"
+{
+ : :: ""
+} 3.238607e-4
Directory: /EOS_Omni/src/
=========================
File [modified]: EOS_Omni_Module.F90
Delta lines: +11 -4
===================================================================
--- EOS_Omni/src/EOS_Omni_Module.F90 2010-07-02 21:41:36 UTC (rev 2)
+++ EOS_Omni/src/EOS_Omni_Module.F90 2010-07-02 21:41:51 UTC (rev 3)
@@ -2,15 +2,19 @@
implicit none
- real*8,parameter :: rho_gf = 1.61930347d-18
- real*8,parameter :: press_gf = 1.80171810d-39
+ real*8,parameter :: rho_gf = 1.61620075314614d-18
+! real*8,parameter :: rho_gf = 1.61930347d-18
+ real*8,parameter :: press_gf = 1.7982953469278d-39
+! real*8,parameter :: press_gf = 1.80171810d-39
real*8,parameter :: eps_gf = 1.11265006d-21
real*8,parameter :: time_gf = 2.03001708d+05
real*8,parameter :: mass_gf = 5.02765209d-34
real*8,parameter :: length_gf = 6.77140812d-06
- real*8,parameter :: inv_rho_gf = 6.17549470205236d17
- real*8,parameter :: inv_press_gf = 5.55025783445257d38
+ real*8,parameter :: inv_rho_gf = 6.18735016707159d17
+! real*8,parameter :: inv_rho_gf = 6.17549470205236d17
+ real*8,parameter :: inv_press_gf = 5.56082181777535d38
+! real*8,parameter :: inv_press_gf = 5.55025783445257d38
real*8,parameter :: inv_eps_gf = 8.98755175549085d20
real*8,parameter :: inv_time_gf = 4.92606692747629d-6
real*8,parameter :: inv_mass_gf = 1.98899999860571d33
@@ -19,4 +23,7 @@
real*8 :: poly_k_cgs = 0.0d0
real*8 :: gl_k_cgs = 0.0d0
+ real*8 :: hybrid_k1_cgs = 0.0d0
+ real*8 :: hybrid_k2_cgs = 0.0d0
+
end module EOS_Omni_Module
File [modified]: EOS_Omni_SingleVarCalls.F90
Delta lines: +106 -13
===================================================================
--- EOS_Omni/src/EOS_Omni_SingleVarCalls.F90 2010-07-02 21:41:36 UTC (rev 2)
+++ EOS_Omni/src/EOS_Omni_SingleVarCalls.F90 2010-07-02 21:41:51 UTC (rev 3)
@@ -24,8 +24,11 @@
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 @@
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 @@
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,7 +144,34 @@
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!"
call CCTK_WARN(0,warnstring)
@@ -138,9 +194,8 @@
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 @@
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 @@
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 @@
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 @@
! 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 @@
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)
File [modified]: EOS_Omni_Startup.F90
Delta lines: +7 -2
===================================================================
--- EOS_Omni/src/EOS_Omni_Startup.F90 2010-07-02 21:41:36 UTC (rev 2)
+++ EOS_Omni/src/EOS_Omni_Startup.F90 2010-07-02 21:41:51 UTC (rev 3)
@@ -11,8 +11,13 @@
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_ARGUMENTS
- poly_k_cgs = poly_k * rho_gf**poly_gamma / press_gf
+ poly_k_cgs = poly_k * rho_gf**poly_gamma_ini / press_gf
- gl_k_cgs = gl_k * rho_gf**gl_gamma / press_gf
+ gl_k_cgs = gl_k * rho_gf**poly_gamma_ini / press_gf
+ hybrid_k1_cgs = hybrid_k1 * rho_gf**poly_gamma_ini / press_gf
+
+ hybrid_k2_cgs = hybrid_k1 * (hybrid_rho_nuc * inv_rho_gf)**(hybrid_gamma1-hybrid_gamma2)
+
+
end subroutine EOS_Omni_Startup
More information about the Commits
mailing list