add EOS Omni support

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@157 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 47 insertions, 3 deletions

diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index 45dfce0..7620ff1 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -253,10 +253,27 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
   CCTK_INT :: type_bits, atmosphere, not_atmosphere
   CCTK_INT :: eos_handle
 
+#if !USE_EOS_OMNI
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+#endif
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+  integer :: n = 1
+  integer :: keytemp = 0
+  integer :: anyerr = 0
+  integer :: keyerr(1) = 0
+  real*8  :: xpress = 0.0d0
+  real*8  :: xeps = 0.0d0
+  real*8  :: xtemp = 0.0d0
+  real*8  :: xye = 0.0d0
+! end EOS Omni vars
+#endif
+
   eos_handle = GRHydro_polytrope_handle
-  if (use_eosgeneral == 0) then
-    eos_handle = -1
-  endif
   
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
@@ -273,6 +290,15 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
           velx(i,j,k) = 0.0d0
           vely(i,j,k) = 0.0d0
           velz(i,j,k) = 0.0d0
+#if USE_EOS_OMNI
+          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
+          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                 rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+#else
+          press(i,j,k) = EOS_Pressure(eos_handle, GRHydro_rho_min, eps(i,j,k))
+          eps(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, GRHydro_rho_min, press(i,j,k))
+#endif
           det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \
                gyy(i,j,k), gyz(i,j,k), gzz(i,j,k))
           call prim2conpolytype(eos_handle, &
@@ -288,6 +314,15 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
             velx_p(i,j,k) = 0.0d0
             vely_p(i,j,k) = 0.0d0
             velz_p(i,j,k) = 0.0d0
+#if USE_EOS_OMNI
+          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                 rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr)
+          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                 rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr)
+#else
+            press_p(i,j,k) = EOS_Pressure(eos_handle, GRHydro_rho_min, eps_p(i,j,k))
+            eps_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, GRHydro_rho_min, press_p(i,j,k))
+#endif
             det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \
                   gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k))
             call prim2conpolytype(eos_handle, &
@@ -304,6 +339,15 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
             velx_p_p(i,j,k) = 0.0d0
             vely_p_p(i,j,k) = 0.0d0
             velz_p_p(i,j,k) = 0.0d0
+#if USE_EOS_OMNI
+          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                 rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr)
+          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                 rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr)
+#else
+            press_p_p(i,j,k) = EOS_Pressure(eos_handle, GRHydro_rho_min, eps_p_p(i,j,k))
+            eps_p_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, GRHydro_rho_min, press_p_p(i,j,k))
+#endif
             det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \
                   gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k))
             call prim2conpolytype(eos_handle, &