RIT MHD dev:

Update initial data routines to properly call out pressures and other
hydro quantities using EOS_Omni rather than now-defunct interfaces.

Add C2P_polytype interface.

Un-hardwire fixed values of gamma for all MHD Con2Prim routines.




git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@117 ac85fae7-cede-4708-beff-ae01c7fa1c26

15 files changed, 363 insertions, 225 deletions

diff --git a/interface.ccl.omni b/interface.ccl.omni
deleted file mode 100644
index e8fe1ae..0000000
--- a/interface.ccl.omni
+++ /dev/null
@@ -1,84 +0,0 @@
-# Interface definition for thorn GRHydro_Init_Data
-# $Header$
-
-implements: GRHydro_init_data
-inherits: GRHydro grid
-
-#USES INCLUDE: EOS_Base.inc
-USES INCLUDE: SpaceMask.h
-#USES INCLUDE: EOS_Base.h
-
-
-SUBROUTINE SpatialDet(CCTK_REAL IN gxx, CCTK_REAL IN gxy, \
-                      CCTK_REAL IN gxz, CCTK_REAL IN gyy, \
-		      CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
-		      CCTK_REAL OUT det)
-
-
-SUBROUTINE Prim2ConPoly(CCTK_INT IN handle, \
-                        CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \
-			CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
-			CCTK_REAL IN det, CCTK_REAL OUT dens, \
-			CCTK_REAL OUT sx, CCTK_REAL OUT sy, \
-			CCTK_REAL OUT sz, CCTK_REAL OUT tau, \
-			CCTK_REAL IN rho, CCTK_REAL IN velx, \
-			CCTK_REAL IN vely, \
-			CCTK_REAL IN velz, CCTK_REAL OUT epsilon, \
-			CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz)
-
-
-SUBROUTINE Prim2ConGen(CCTK_INT IN handle, \
-                          CCTK_REAL IN gxx, CCTK_REAL IN gxy, \
-                          CCTK_REAL IN gxz, CCTK_REAL IN gyy, \
-                          CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
-                          CCTK_REAL IN det, CCTK_REAL OUT dens, \
-                          CCTK_REAL OUT sx, CCTK_REAL OUT sy, \
-                          CCTK_REAL OUT sz, CCTK_REAL OUT tau, \
-                          CCTK_REAL IN rho, CCTK_REAL IN velx, \
-                          CCTK_REAL IN vely, \
-                          CCTK_REAL IN velz, CCTK_REAL IN epsilon, \
-                          CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz)
-
-SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \
-                        CCTK_REAL OUT sx, CCTK_REAL OUT sy, \
-			CCTK_REAL OUT sz, CCTK_REAL OUT tau, \
-			CCTK_REAL OUT rho, CCTK_REAL OUT velx, \
-			CCTK_REAL OUT vely, CCTK_REAL OUT velz, \
-			CCTK_REAL OUT epsilon, CCTK_REAL OUT press, \
-			CCTK_REAL OUT w_lorentz, CCTK_REAL IN uxx, \
-			CCTK_REAL IN uxy, CCTK_REAL IN uxz, CCTK_REAL IN uyy, \
-			CCTK_REAL IN uyz, CCTK_REAL IN uzz, CCTK_REAL IN det, \
-			CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \
-			CCTK_REAL IN r, CCTK_REAL IN rho_min, \
-			CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed)
-
-
-USES FUNCTION SpatialDet
-USES FUNCTION Prim2ConPoly
-USES FUNCTION Prim2ConGen
-USES FUNCTION Con2PrimPoly
-
-protected:
-
-CCTK_REAL simple_wave_grid_functions TYPE=GF TIMELEVELS=1 tags='checkpoint="no"'
-{
-  simple_tmp
-  c_s
-} "1D arrays for the simple-wave routine"
-
-CCTK_REAL simple_wave_scalars TYPE=scalar
-{
-     simple_rho_0
-     simple_eps_0
-} "values at v=0"
-
-CCTK_REAL simple_wave_output TYPE=GF TIMELEVELS=1 tags='checkpoint="no"'
-{
-  simple_rho
-  simple_eps
-#  simple_entropy
-} "output variables for the simple-wave routine"
-
-private:
-
-CCTK_INT GRHydro_init_data_reflevel type = SCALAR tags='checkpoint="no"' "Refinement level GRHydro is working on right now"
diff --git a/param.ccl b/param.ccl
index bc000b0..01fbe7f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -22,6 +22,7 @@ EXTENDS KEYWORD initial_data ""
   "con2primtest"	:: "Testing the con -> prim conversion"
   "con2prim2con_test"	:: "Testing the con -> prim -> con conversion"
   "prim2con2prim_test"	:: "Testing the prim -> con -> prim conversion"
+  "prim2con2prim_polytype_test"	:: "Testing the prim -> con -> prim conversion - polytype version"
   "reconstruction_test"	:: "Testing reconstruction"
 }
 
@@ -41,6 +42,11 @@ KEYWORD shock_case "Simple, Sod's problem or other?"
   "Simple"	:: "GRAstro_Hydro test case"
   "Sod"		:: "Sod's problem"
   "Blast"	:: "Strong blast wave"
+  "Balsara1"    :: "Balsara Test #1"
+  "Balsara2"    :: "Balsara Test #2"
+  "Balsara3"    :: "Balsara Test #3"
+  "Balsara4"    :: "Balsara Test #4"
+  "Balsara5"    :: "Balsara Test #5"
 } "Sod"
 
 REAL shock_xpos "Position of shock plane: x"
@@ -106,6 +112,7 @@ REAL Bz_init "Initial B-field in the z-dir"
 
 shares:GRHydro
 
+USES real GRHydro_eos_rf_prec
 USES real GRHydro_rho_central
 USES real GRHydro_eps_min ""
 USES real GRHydro_perc_ptol ""
diff --git a/schedule.ccl b/schedule.ccl
index 26bdee4..a35dfdb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -77,6 +77,22 @@ if (CCTK_Equals(initial_data,"prim2con2prim_test")) {
   }
 }
 
+if (CCTK_Equals(initial_data,"prim2con2prim_polytype_test")) {
+  STORAGE:GRHydro_init_data_reflevel
+  schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE p2c2p_call
+  {
+    LANG: Fortran
+  } "Calculate current refinement level"
+
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
+  {
+    schedule p2c2pM_polytype in HydroBase_Initial AS p2c2p_call
+    {
+      LANG: Fortran
+    } "Testing primitive to conservative to primitive - MHD polytype version"
+  }
+}
+
 if (CCTK_Equals(initial_data,"reconstruction_test")) {
   schedule GRHydro_reconstruction_test in HydroBase_Initial
   {
diff --git a/src/CheckParam.c b/src/CheckParam.c
index 31222eb..124db7f 100644
--- a/src/CheckParam.c
+++ b/src/CheckParam.c
@@ -23,6 +23,15 @@ void GRHydro_InitData_CheckParameters(CCTK_ARGUMENTS)
     {
       CCTK_PARAMWARN("That test not yet implemented in MHD!");
     }
-	  
+
+  if(!CCTK_Equals(Bvec_evolution_method,"GRHydro") && 
+     ((CCTK_Equals(shock_case,"Balsara1")) ||
+      (CCTK_Equals(shock_case,"Balsara2")) ||
+      (CCTK_Equals(shock_case,"Balsara3")) ||
+      (CCTK_Equals(shock_case,"Balsara4")) ||
+      (CCTK_Equals(shock_case,"Balsara5"))))
+    {
+      CCTK_PARAMWARN("That test requires MHD!  Set Bvec_evolution_method=GRHYDRO!");
+    }
 }
 
diff --git a/src/GRHydro_C2P2C.F90 b/src/GRHydro_C2P2C.F90
index 1ce1347..a1cc36d 100644
--- a/src/GRHydro_C2P2C.F90
+++ b/src/GRHydro_C2P2C.F90
@@ -34,13 +34,6 @@ subroutine c2p2c(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   integer didit,i,j,k,nx,ny,nz
   CCTK_REAL det
   CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
@@ -52,7 +45,6 @@ subroutine c2p2c(CCTK_ARGUMENTS)
   CCTK_INT C2P_failed
   logical epsnegative
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -63,7 +55,6 @@ subroutine c2p2c(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
 
@@ -105,16 +96,13 @@ subroutine c2p2c(CCTK_ARGUMENTS)
   epsnegative = .false.
   
   GRHydro_rho_min = 1.e-10
-#if USE_EOS_OMNI
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
 
-  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  pmin   = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
-#endif
+
   C2P_failed = 0
 
   write(*,*) 'C2P2C test: initial values.'
diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90
index 12c1ed4..3bd61b6 100644
--- a/src/GRHydro_C2P2CM.F90
+++ b/src/GRHydro_C2P2CM.F90
@@ -34,13 +34,6 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   integer didit,i,j,k,nx,ny,nz
   CCTK_REAL det
   CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
@@ -49,11 +42,10 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send
   CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send
   CCTK_REAL bvcx_send, bvcy_send, bvcz_send, b2_send
-  CCTK_REAL pmin, epsmin
+  CCTK_REAL pmin, epsmin, local_gam
   CCTK_INT C2P_failed
   logical epsnegative
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -64,7 +56,6 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
 
@@ -110,16 +101,18 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   epsnegative = .false.
   
   GRHydro_rho_min = 1.e-10
-#if USE_EOS_OMNI
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
 
-  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  pmin   = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
-#endif
+
+  local_gam=0.0
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       1.d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam = local_gam+1.0
+
   C2P_failed = 0
 
   write(*,*) 'C2P2CM test: initial values.'
@@ -136,7 +129,7 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   write(*,*) '   Bvecz  : ',bvcz_send
   
   write(*,*) 'C2P2CM test: getting the associated primitive variables.'
-  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,dens_send,sx_send,sy_send,sz_send, &
+  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,dens_send,sx_send,sy_send,sz_send, &
        tau_send,rho_send,velx_send,vely_send,velz_send, &
        eps_send,press_send,w_lorentz_send, &
        gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,&
diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index d6950ab..6b9ba34 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -52,14 +52,6 @@ subroutine GRHydro_con2primtest(CCTK_ARGUMENTS)
   CCTK_INT C2P_failed
   logical epsnegative
 
-#ifndef 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
@@ -70,7 +62,6 @@ subroutine GRHydro_con2primtest(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call CCTK_WARN(1,"For this test, remember to use a polytropic EoS and to set eos_gamma = 2.0 and eos_k = 100.0")
 
@@ -107,16 +98,12 @@ subroutine GRHydro_con2primtest(CCTK_ARGUMENTS)
 
   epsnegative = .false.
 
-#if USE_EOS_OMNI
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
-       GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr)
 
-  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  pmin   = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
-#endif
+
   C2P_failed = 0
 
   write(*,*) 'Con2Prim test: converting to primitive variables'
diff --git a/src/GRHydro_P2C2P.F90 b/src/GRHydro_P2C2P.F90
index 10258c3..ba13d12 100644
--- a/src/GRHydro_P2C2P.F90
+++ b/src/GRHydro_P2C2P.F90
@@ -34,13 +34,6 @@ subroutine p2c2p(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   integer didit,i,j,k,nx,ny,nz
   CCTK_REAL det
   CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
@@ -52,7 +45,6 @@ subroutine p2c2p(CCTK_ARGUMENTS)
   CCTK_INT C2P_failed
   logical epsnegative
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -63,7 +55,6 @@ subroutine p2c2p(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
 
@@ -103,19 +94,15 @@ subroutine p2c2p(CCTK_ARGUMENTS)
   epsnegative = .false.
   
   GRHydro_rho_min = 1.e-10
-#if USE_EOS_OMNI
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
-       rho_send,1.0d0,xtemp,eps_send,press_send,keyerr,anyerr)
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       rho_send,eps_send,xtemp,xye,press_send,keyerr,anyerr)
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
 
-  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  press_send   = EOS_Pressure(GRHydro_eos_handle, rho_send, eps_send)
-  pmin   = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
-#endif
+
   C2P_failed = 0
 
   write(*,*) 'P2C2P test: the primitive variables are'
diff --git a/src/GRHydro_P2C2PM.F90 b/src/GRHydro_P2C2PM.F90
index c68b735..1abf1ca 100644
--- a/src/GRHydro_P2C2PM.F90
+++ b/src/GRHydro_P2C2PM.F90
@@ -34,13 +34,6 @@ subroutine p2c2pm(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   integer didit,i,j,k,nx,ny,nz
   CCTK_REAL det
   CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
@@ -49,11 +42,10 @@ subroutine p2c2pm(CCTK_ARGUMENTS)
   CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send
   CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send
   CCTK_REAL bvcx_send,bvcy_send,bvcz_send,b2_send
-  CCTK_REAL pmin, epsmin
+  CCTK_REAL pmin, epsmin, local_gam
   CCTK_INT C2P_failed
   logical epsnegative
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -64,7 +56,6 @@ subroutine p2c2pm(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
 
@@ -108,19 +99,19 @@ subroutine p2c2pm(CCTK_ARGUMENTS)
   epsnegative = .false.
   
   GRHydro_rho_min = 1.e-10
-#if USE_EOS_OMNI
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
-       rho_send,1.0d0,xtemp,eps_send,press_send,keyerr,anyerr)
-  call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       rho_send,eps_send,xtemp,xye,press_send,keyerr,anyerr)
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
 
-  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  press_send   = EOS_Pressure(GRHydro_eos_handle, rho_send, eps_send)
-  pmin   = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
-#endif
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       1.0d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam=local_gam+1.0
+
   C2P_failed = 0
 
   write(*,*) 'P2C2PM test: the primitive variables are'
@@ -155,7 +146,8 @@ subroutine p2c2pm(CCTK_ARGUMENTS)
   write(*,*) '   Bvecz  : ',bvcz_send
   
   write(*,*) 'P2C2PM test: getting the associated primitive variables.'
-  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,dens_send,sx_send,sy_send,sz_send, &
+  
+  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,dens_send,sx_send,sy_send,sz_send, &
        tau_send,rho_send,velx_send,vely_send,velz_send, &
        eps_send,press_send,w_lorentz_send, &
        gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,&
diff --git a/src/GRHydro_P2C2PM_polytype.F90 b/src/GRHydro_P2C2PM_polytype.F90
new file mode 100644
index 0000000..a8d0e77
--- /dev/null
+++ b/src/GRHydro_P2C2PM_polytype.F90
@@ -0,0 +1,181 @@
+ /*@@
+   @file      GRHydro_P2C2PM_polytype.F90
+   @date      Sep 25, 2010
+   @author    Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+   @desc 
+   A test of the conservative <--> primitive variable exchange
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+ /*@@
+   @routine    p2c2pm
+   @date       Sep 25, 2010
+   @author     Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+   @desc 
+   Testing the conservative <--> primitive variable transformations.
+   The values before and after should match.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine p2c2pm_polytype(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  integer didit,i,j,k,nx,ny,nz
+  CCTK_REAL det, local_K
+  CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
+  CCTK_REAL gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send
+  CCTK_REAL dens_send,sx_send,sy_send,sz_send,tau_send
+  CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send
+  CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send
+  CCTK_REAL bvcx_send,bvcy_send,bvcz_send,b2_send
+  CCTK_REAL pmin, epsmin, local_gam, sc_send
+  CCTK_INT C2P_failed
+  logical epsnegative
+
+! 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
+
+  call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  x_send = 0.0d0
+  y_send = 0.0d0
+  z_send = 0.0d0
+  r_send = 0.0d0
+  
+  gxx_send = 1.0d0
+  gyy_send = 1.0d0 
+  gzz_send = 1.0d0
+  gxy_send = 0.0d0
+  gxz_send = 0.0d0
+  gyz_send = 0.0d0
+  
+  det = 1.0d0
+  
+  uxx = 1.0d0
+  uyy = 1.0d0 
+  uzz = 1.0d0
+  uxy = 0.0d0
+  uxz = 0.0d0
+  uyz = 0.0d0
+  
+  rho_send = 1.29047362d0
+  velx_send   = 0.166666658d0
+  vely_send   = 0.166666658d0
+  velz_send   = 0.166666658d0
+!!$  eps_send  = 0.484123939d0
+  
+  bvcx_send   = Bx_init
+  bvcy_send   = By_init
+  bvcz_send   = Bz_init
+
+  w_lorentz_send = 1.d0/sqrt(1.0d0-velx_send*velx_send-vely_send*vely_send-velz_send*velz_send)
+
+  epsnegative = .false.
+  
+  GRHydro_rho_min = 1.e-10
+
+  call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       rho_send,1.0d0,xtemp,xye,press_send,keyerr,anyerr)
+  call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+  call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       rho_send,xeps,xtemp,xye,press_send,eps_send,keyerr,anyerr)
+  call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+
+  call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       1.d0,1.d0,xtemp,xye,xpress,keyerr,anyerr)
+  local_K = xpress
+
+  local_gam=log(press_send/local_K)/log(rho_send)
+
+  C2P_failed = 0
+
+  write(*,*) 'P2C2PM test: the primitive variables are'
+  write(*,*) '   primitive variables: '
+  write(*,*) '   rho   : ',rho_send
+  write(*,*) '   velx  : ',velx_send
+  write(*,*) '   vely  : ',vely_send
+  write(*,*) '   velz  : ',velz_send
+  write(*,*) '   press : ',press_send
+  write(*,*) '   eps   : ',eps_send
+  write(*,*) '   W     : ',w_lorentz_send
+  write(*,*) '   Bvecx  : ',bvcx_send
+  write(*,*) '   Bvecy  : ',bvcy_send
+  write(*,*) '   Bvecz  : ',bvcz_send
+  
+  write(*,*) 'P2C2PM test: converting back to conserved variables.'
+  call prim2conpolytypeM(GRHydro_polytrope_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+       dens_send, sx_send, sy_send, sz_send, tau_send, bvcx_send, bvcy_send, bvcz_send, rho_send, &
+       velx_send, vely_send, velz_send, eps_send, press_send, w_lorentz_send) 
+
+  write(*,*) 'P2C2PM test: initial values.'
+  write(*,*) '   conservative variables: '
+  write(*,*) '   dens: ',dens_send
+  write(*,*) '   sx  : ',sx_send
+  write(*,*) '   sy  : ',sy_send
+  write(*,*) '   sz  : ',sz_send
+  write(*,*) '   tau : ',tau_send
+  write(*,*) '   eps : ',eps_send
+  write(*,*) '   W   : ',w_lorentz_send
+  write(*,*) '   Bvecx  : ',bvcx_send
+  write(*,*) '   Bvecy  : ',bvcy_send
+  write(*,*) '   Bvecz  : ',bvcz_send
+
+  write(6,*)'local_K:',local_K
+  sc_send=local_K*dens_send
+
+  write(*,*) 'P2C2PM test: getting the associated primitive variables.'
+  call GRHydro_Con2PrimM_Polytype_pt(GRHydro_polytrope_handle,local_gam,dens_send,sx_send,sy_send,sz_send, &
+       sc_send,rho_send,velx_send,vely_send,velz_send, &
+       eps_send,press_send,w_lorentz_send, &
+       gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,&
+       uxx,uxy,uxz,uyy,uyz,uzz,det,&
+       bvcx_send,bvcy_send,bvcz_send,b2_send,&
+       epsnegative,C2P_failed)
+  
+  write(*,*) 'P2C2PM test: the primitive variables are'
+  write(*,*) '   primitive variables: '
+  write(*,*) '   rho   : ',rho_send
+  write(*,*) '   velx  : ',velx_send
+  write(*,*) '   vely  : ',vely_send
+  write(*,*) '   velz  : ',velz_send
+  write(*,*) '   press : ',press_send
+  write(*,*) '   eps   : ',eps_send
+  write(*,*) '   W     : ',w_lorentz_send
+  write(*,*) '   Bvecx  : ',bvcx_send
+  write(*,*) '   Bvecy  : ',bvcy_send
+  write(*,*) '   Bvecz  : ',bvcz_send
+  
+  STOP
+
+  return
+
+end subroutine p2c2pm_polytype
diff --git a/src/GRHydro_ReadConformalData.F90 b/src/GRHydro_ReadConformalData.F90
index c492d3c..665b91e 100644
--- a/src/GRHydro_ReadConformalData.F90
+++ b/src/GRHydro_ReadConformalData.F90
@@ -43,18 +43,9 @@ subroutine GRHydro_ReadConformalData(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
-!#include "EOS_Base.h"
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-
-  
   CCTK_INT :: i,j,k,handle,ierr
-!  CCTK_REAL :: eos_k, eos_gamma
   CCTK_REAL :: rho_min, det
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: poly_eoskey = 0
@@ -68,7 +59,6 @@ subroutine GRHydro_ReadConformalData(CCTK_ARGUMENTS)
   real*8  :: xye(1) = 0.0d0
   poly_eoskey = EOS_Omni_GetHandle("2D_Polytrope")
 ! end EOS Omni vars
-#endif
 
   ! only gxx has been read in; copy it into gyy and gzz as well
   gyy = gxx
@@ -155,14 +145,11 @@ subroutine GRHydro_ReadConformalData(CCTK_ARGUMENTS)
   do i=1,cctk_lsh(1)
     do j=1,cctk_lsh(2)
       do k=1,cctk_lsh(3)
-#ifdef USE_EOS_OMNI
+
         call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,rf_precision,&
              n,rho(i,j,k),xeps,xtemp,xye,&
              press(i,j,k),eps(i,j,k),keyerr,anyerr)
 
-#else
-        eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_eos_handle,rho(i,j,k),press(i,j,k))
-#endif
       end do
     end do
   end do
diff --git a/src/GRHydro_ShockTube.F90 b/src/GRHydro_ShockTube.F90
index 80d8968..baf891d 100644
--- a/src/GRHydro_ShockTube.F90
+++ b/src/GRHydro_ShockTube.F90
@@ -11,6 +11,7 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 #include "cctk_Functions.h"
+#include "GRHydro_Macros.h"
 
 #define velx(i,j,k) vel(i,j,k,1)
 #define vely(i,j,k) vel(i,j,k,2)
@@ -124,8 +125,7 @@ subroutine GRHydro_shocktube(CCTK_ARGUMENTS)
           eps(i,j,k) = epsr
         end if
 
-        call SpatialDet(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
-             gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det)
+        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))
 
         if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then
           call Prim2ConPoly(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90
index 1169fef..96b0e7c 100644
--- a/src/GRHydro_ShockTubeM.F90
+++ b/src/GRHydro_ShockTubeM.F90
@@ -54,6 +54,13 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
        velzl, velzr, epsl, epsr
   CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr
   
+  bvcxl = Bx_init
+  bvcyl = By_init
+  bvczl = Bz_init
+  bvcxr = Bx_init
+  bvcyr = By_init
+  bvczr = Bz_init
+
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
   nz = cctk_lsh(3)
@@ -110,17 +117,95 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           velzr = 0.d0
           epsl = 1500.d0
           epsr = 1.5d-2
-        else
+       else if (CCTK_EQUALS(shock_case,"Balsara1")) then
+          rhol = 1.d0
+          rhor = 0.125d0
+          velxl = 0.d0
+          velxr = 0.d0
+          velyl = 0.d0
+          velyr = 0.d0
+          velzl = 0.d0
+          velzr = 0.d0
+          bvcxl=0.5d0
+          bvcxr=0.5d0
+          bvcyl=1.d0
+          bvcyr=-1.d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 1.0/rhol
+          epsr = 0.1/rhor
+       else if (CCTK_EQUALS(shock_case,"Balsara2")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          velxl = 0.d0
+          velxr = 0.d0
+          velyl = 0.d0
+          velyr = 0.d0
+          velzl = 0.d0
+          velzr = 0.d0
+          bvcxl=5.0d0
+          bvcxr=5.0d0
+          bvcyl=6.d0
+          bvcyr=0.7d0
+          bvczl=6.d0
+          bvczr=0.7d0
+          epsl = 1.5d0*30.0d0/rhol
+          epsr = 1.5d0*1.0d0/rhor
+       else if (CCTK_EQUALS(shock_case,"Balsara3")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          velxl = 0.d0
+          velxr = 0.d0
+          velyl = 0.d0
+          velyr = 0.d0
+          velzl = 0.d0
+          velzr = 0.d0
+          bvcxl=10.0d0
+          bvcxr=10.0d0
+          bvcyl=7.d0
+          bvcyr=0.7d0
+          bvczl=7.d0
+          bvczr=0.7d0
+          epsl = 1.5d0*1000.0d0/rhol
+          epsr = 1.5d0*0.1d0/rhor
+       else if (CCTK_EQUALS(shock_case,"Balsara4")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          velxl = 0.999d0
+          velxr = -0.999d0
+          velyl = 0.d0
+          velyr = 0.d0
+          velzl = 0.d0
+          velzr = 0.d0
+          bvcxl=10.0d0
+          bvcxr=10.0d0
+          bvcyl=7.d0
+          bvcyr=-7.d0
+          bvczl=7.d0
+          bvczr=-7.d0
+          epsl = 1.5d0*0.1d0/rhol
+          epsr = 1.5d0*0.1d0/rhor
+       else if (CCTK_EQUALS(shock_case,"Balsara5")) then
+          rhol = 1.08d0
+          rhor = 1.d0
+          velxl = 0.4d0
+          velxr = -0.45d0
+          velyl = 0.3d0
+          velyr = -0.2d0
+          velzl = 0.2d0
+          velzr = 0.2d0
+          bvcxl=2.0d0
+          bvcxr=2.0d0
+          bvcyl=0.3d0
+          bvcyr=-0.7d0
+          bvczl=0.3d0
+          bvczr=0.5d0
+          epsl = 1.5d0*0.95d0/rhol
+          epsr = 1.5d0*1.0d0/rhor
+       else
           call CCTK_WARN(0,"Shock case not recognized")
         end if
 
-        bvcxl = Bx_init
-        bvcyl = By_init
-        bvczl = Bz_init
-        bvcxr = Bx_init
-        bvcyr = By_init
-        bvczr = Bz_init
-
         if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& 
              ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then
           rho(i,j,k) = rhol
diff --git a/src/GRHydro_SimpleWave.F90 b/src/GRHydro_SimpleWave.F90
index 2991ff4..c13916b 100644
--- a/src/GRHydro_SimpleWave.F90
+++ b/src/GRHydro_SimpleWave.F90
@@ -46,17 +46,9 @@ subroutine GRHydro_SimpleWave(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
   
-#ifndef USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   CCTK_INT :: i, j, k, nx, ny, nz
   CCTK_REAL :: dr, k1, k2, k3, k4, in_data, old_data, source_data, new_data, c_0, det, pi
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -68,7 +60,6 @@ subroutine GRHydro_SimpleWave(CCTK_ARGUMENTS)
   real*8  :: xye(1) = 0.0d0
   real*8  :: rf_precision = 1.0d-10
 ! end EOS Omni vars
-#endif
   
   call CCTK_INFO("Setting initial data for a simple wave as Anile Miller Motta")
 
@@ -143,21 +134,19 @@ subroutine GRHydro_SimpleWave(CCTK_ARGUMENTS)
     
     ! atmosphere
 
-    if ( (rho(i,1,1) < GRHydro_rho_min).OR.(velx(i,1,1) < 0) ) then
-      rho(i,1,1) = rho_abs_min
-!      rho(i,1,1) = 1.0 !the value of rho_min for the initial data
-      eps(i,1,1) = rho_abs_min**(1.d0/3.d0)
-      velx(i,1,1) = 0.d0
-      w_lorentz(i,1,1) = 1.d0
-#if USE_EOS_OMNI
-      xeps = 1.0d0
-  call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,rf_precision,n,&
-       rho(i,1,1),xeps,xtemp,xye,press(i,1,1),keyerr,anyerr)
-#else
-      press(i,1,1) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,1,1), 1.0d0)  
-#endif
-      ! polytrope only (initial data)
-    end if
+     if ( (rho(i,1,1) < GRHydro_rho_min).OR.(velx(i,1,1) < 0) ) then
+        rho(i,1,1) = rho_abs_min
+        !      rho(i,1,1) = 1.0 !the value of rho_min for the initial data
+        eps(i,1,1) = rho_abs_min**(1.d0/3.d0)
+        velx(i,1,1) = 0.d0
+        w_lorentz(i,1,1) = 1.d0
+        
+        xeps = 1.0d0
+        call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,rf_precision,n,&
+             rho(i,1,1),xeps,xtemp,xye,press(i,1,1),keyerr,anyerr)
+        
+        ! polytrope only (initial data)
+     end if
 
 !    write(*,*) 'p',i, x(i,1,1), rho(i,1,1)**(4.d0/3.d0)/press(i,1,1)    
     
diff --git a/src/make.code.defn b/src/make.code.defn
index dc78381..09a09f2 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -13,7 +13,8 @@ SRCS = 	GRHydro_C2P2C.F90 \
 	GRHydro_SimpleWave.F90 \
    CheckParam.c\
    GRHydro_P2C2P.F90 \
-   GRHydro_P2C2PM.F90
+   GRHydro_P2C2PM.F90 \
+   GRHydro_P2C2PM_polytype.F90
 
 # Subdirectories containing source files
 SUBDIRS =