make new EOS call that returns press and cs2 [and eps if needed]

From: Christian D. Ott <cott@tapir.caltech.edu>

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

2 files changed, 192 insertions, 8 deletions

diff --git a/interface.ccl b/interface.ccl
index 44a2f11..5b4ef5b 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -22,6 +22,21 @@ void FUNCTION EOS_Omni_press(CCTK_INT IN eoskey,         \
 
 PROVIDES FUNCTION EOS_Omni_press WITH EOS_Omni_EOS_Press LANGUAGE Fortran
 
+void FUNCTION EOS_Omni_press_cs2(CCTK_INT IN eoskey,         \
+			     CCTK_INT IN havetemp,       \
+			     CCTK_REAL IN rf_precision,  \
+                             CCTK_INT IN npoints,        \
+			     CCTK_REAL IN ARRAY rho,     \
+			     CCTK_REAL INOUT ARRAY eps,  \
+			     CCTK_REAL INOUT ARRAY temp, \
+			     CCTK_REAL IN ARRAY ye,      \
+			     CCTK_REAL OUT ARRAY press,  \
+			     CCTK_REAL OUT ARRAY cs2,    \
+			     CCTK_INT OUT ARRAY keyerr,  \
+                             CCTK_INT OUT anyerr)
+
+PROVIDES FUNCTION EOS_Omni_press_cs2 WITH EOS_Omni_EOS_Press_cs2 LANGUAGE Fortran
+
 void FUNCTION EOS_Omni_pressOMP(CCTK_INT IN eoskey,         \
 			     CCTK_INT IN havetemp,       \
 			     CCTK_REAL IN rf_precision,  \
diff --git a/src/EOS_Omni_MultiVarCalls.F90 b/src/EOS_Omni_MultiVarCalls.F90
index 2b3ea1b..379f8db 100644
--- a/src/EOS_Omni_MultiVarCalls.F90
+++ b/src/EOS_Omni_MultiVarCalls.F90
@@ -101,7 +101,6 @@ subroutine EOS_Omni_EOS_full(eoskey,keytemp,rf_precision,npoints,&
   CCTK_REAL, intent(out)   :: muhat(npoints)
 
   ! local vars
-  integer          :: i
   character(256)   :: warnstring
 
   if(eoskey.ne.4) then
@@ -157,10 +156,6 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: hybrid_local_gamma
   real*8           :: hybrid_local_k_cgs
   real*8           :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2
-  ! temporary vars for nuc_eos
-  real*8           :: xrho,xye,xtemp,xenr,xent
-  real*8           :: xprs,xmunu,xcs2
-  real*8           :: xdedt,xdpderho,xdpdrhoe
 
   anyerr    = 0
   keyerr(:) = 0
@@ -264,9 +259,7 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: hybrid_local_k_cgs
   real*8           :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2
   ! temporary vars for nuc_eos
-  real*8           :: xrho,xye,xtemp,xenr,xent
-  real*8           :: xprs,xmunu,xcs2
-  real*8           :: xdedt,xdpderho,xdpdrhoe
+  real*8           :: xdpderho,xdpdrhoe
 
 
   anyerr    = 0
@@ -342,3 +335,179 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
    end select
 
  end subroutine EOS_Omni_EOS_DEpsByDRho_DEpsByDPress
+
+
+subroutine EOS_Omni_EOS_Press_cs2(eoskey,keytemp,rf_precision,npoints,&
+                              rho,eps,temp,ye,press,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)   :: press(npoints)
+  CCTK_REAL, intent(out)   :: cs2(npoints)
+
+  ! local vars
+  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
+  ! temporary vars for nuc_eos
+  real*8           :: xrho
+  real*8           :: 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, 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) = press_gf * poly_k_cgs * &
+                   (rho(i)*inv_rho_gf)**(poly_gamma) / &
+                   (poly_gamma - 1.0d0) / rho(i)
+           enddo
+        endif
+        do i=1,npoints
+           press(i) = press_gf * poly_k_cgs * &
+                (rho(i)*inv_rho_gf)**poly_gamma
+           cs2(i) = poly_gamma * press(i) / rho(i) / &
+                (1 + eps(i) + press(i)/rho(i))
+        enddo
+     case (2)
+        ! gamma-law EOS
+        if(keytemp.eq.1) then
+           do i=1,npoints
+              eps(i) = press_gf * gl_k_cgs * &
+                   (rho(i)*inv_rho_gf)**(gl_gamma) / &
+                   (gl_gamma - 1.0d0) / rho(i)
+           enddo
+        endif
+        do i=1,npoints
+           press(i) = (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 + press(i) * xdpderho / (rho(i)**2)) / &
+                (1.0d0 + eps(i) + press(i)/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
+           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
+
+           cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / &
+                rho(i) / (1.0d0 + eps(i) + press(i)/rho(i))
+        enddo
+
+     case (4)
+        ! nuc eos
+        if(keytemp.eq.1) then
+           call nuc_eos_m_kt1_press_eps_cs2(npoints,rho,temp,ye,&
+                eps,press,cs2,keyerr,anyerr)
+        else if(keytemp.eq.0) then
+           call nuc_eos_m_kt0_press_cs2(npoints,rho,temp,ye,eps,press,&
+                cs2,rf_precision,keyerr,anyerr)
+        else
+           call CCTK_ERROR("This keytemp is not suppported!")
+           STOP
+        endif
+     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
+              press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma
+              eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i) 
+              cs2(i) = coldeos_low_kappa * coldeos_low_gamma * &
+                   rho(i)**(coldeos_low_gamma-1.0d0) / &
+                   (1.0d0 + eps(i) + press(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)
+
+           ! 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
+           press(i) = press_cold + press_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 (6)
+        ! barotropic tabular EOS with gamma law
+        write(warnstring,*) "eoskey ",eoskey," for barotropic EOS not implemented!"
+        call CCTK_ERROR(warnstring)
+        STOP
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_ERROR(warnstring)
+        STOP
+     end select
+
+end subroutine EOS_Omni_EOS_Press_cs2
+