* add hybrid EOS

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/EOS_Omni@3 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

4 files changed, 157 insertions, 21 deletions

diff --git a/param.ccl b/param.ccl
index 6b8acf7..870d04d 100644
--- a/param.ccl
+++ b/param.ccl
@@ -4,7 +4,12 @@
 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
@@ -15,7 +20,7 @@ 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 @@ REAL gl_k "Polytropic constant in c=G=Msun=1 for gamma-law EOS"
 {
  : :: ""
 } 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
diff --git a/src/EOS_Omni_Module.F90 b/src/EOS_Omni_Module.F90
index 702d841..75c6c2d 100644
--- a/src/EOS_Omni_Module.F90
+++ b/src/EOS_Omni_Module.F90
@@ -2,15 +2,19 @@ module EOS_Omni_Module
 
   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 @@ module EOS_Omni_Module
   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
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)
diff --git a/src/EOS_Omni_Startup.F90 b/src/EOS_Omni_Startup.F90
index 7a9a2f7..c2a5db9 100644
--- a/src/EOS_Omni_Startup.F90
+++ b/src/EOS_Omni_Startup.F90
@@ -11,8 +11,13 @@ subroutine EOS_Omni_Startup(CCTK_ARGUMENTS)
   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**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)
 
-  gl_k_cgs   = gl_k * rho_gf**gl_gamma / press_gf
 
 end subroutine EOS_Omni_Startup