Add the magnetized routine versions for prim2con2prim

and vice versa. 

Add the magnetized version for the ShockTube routine.



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

10 files changed, 827 insertions, 7 deletions

diff --git a/interface.ccl b/interface.ccl
index 70912bd..b8f7af7 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -39,6 +39,32 @@ SUBROUTINE Prim2ConGen(CCTK_INT IN handle, \
                           CCTK_REAL IN velz, CCTK_REAL IN epsilon, \
                           CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz)
 
+SUBROUTINE Prim2ConPolyM(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 Bvecx, CCTK_REAL IN Bvecy, \
+			CCTK_REAL IN Bvecz, 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 Prim2ConGenM(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 Bvecx, CCTK_REAL IN Bvecy, \
+			  CCTK_REAL IN Bvecz, 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, \
@@ -56,6 +82,8 @@ SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \
 USES FUNCTION SpatialDet
 USES FUNCTION Prim2ConPoly
 USES FUNCTION Prim2ConGen
+USES FUNCTION Prim2ConPolyM
+USES FUNCTION Prim2ConGenM
 USES FUNCTION Con2PrimPoly
 
 CCTK_INT FUNCTION EOS_Omni_GetHandle(CCTK_STRING IN name)
diff --git a/param.ccl b/param.ccl
index 2655115..bc000b0 100644
--- a/param.ccl
+++ b/param.ccl
@@ -4,6 +4,7 @@
 shares:HydroBase
 
 USES CCTK_INT timelevels
+USES KEYWORD Bvec_evolution_method
 
 EXTENDS KEYWORD initial_hydro ""
 {
@@ -20,6 +21,7 @@ EXTENDS KEYWORD initial_data ""
 #  "shocktube"		:: "Shock tube initial data for GRHydro"
   "con2primtest"	:: "Testing the con -> prim conversion"
   "con2prim2con_test"	:: "Testing the con -> prim -> con conversion"
+  "prim2con2prim_test"	:: "Testing the prim -> con -> prim conversion"
   "reconstruction_test"	:: "Testing reconstruction"
 }
 
@@ -86,6 +88,22 @@ BOOLEAN attenuate_atmosphere "Attenuate the velocity in the atmosphere"
 {
 } "no"
 
+REAL Bx_init "Initial B-field in the x-dir"
+{
+  *:*     :: "Anything"
+} 0.0
+
+REAL By_init "Initial B-field in the y-dir"
+{
+  *:*     :: "Anything"
+} 0.0
+
+REAL Bz_init "Initial B-field in the z-dir"
+{
+  *:*     :: "Anything"
+} 0.0
+
+
 shares:GRHydro
 
 USES real GRHydro_rho_central
diff --git a/schedule.ccl b/schedule.ccl
index 1622776..26bdee4 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -6,10 +6,20 @@ schedule GRHydro_InitData_CheckParameters  AT CCTK_PARAMCHECK
 } "Check parameters"
 
 if (CCTK_Equals(initial_hydro,"shocktube")) {
-  schedule GRHydro_shocktube in HydroBase_Initial
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
   {
-    LANG: Fortran
-  } "Shocktube initial data"
+    schedule GRHydro_shocktubeM in HydroBase_Initial
+    {
+      LANG: Fortran
+    } "Shocktube initial data - MHD version"
+
+  } else {
+    schedule GRHydro_shocktube in HydroBase_Initial
+    {
+      LANG: Fortran
+    } "Shocktube initial data"
+
+  }
 }
 
 if (CCTK_Equals(initial_data,"con2primtest")) {
@@ -27,15 +37,44 @@ if (CCTK_Equals(initial_data,"con2primtest")) {
 
 if (CCTK_Equals(initial_data,"con2prim2con_test")) {
   STORAGE:GRHydro_init_data_reflevel
-  schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE c2p2c
+  schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE c2p2c_call
   {
     LANG: Fortran
   } "Calculate current refinement level"
 
-  schedule c2p2c in HydroBase_Initial
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
+  {
+    schedule c2p2cM in HydroBase_Initial AS c2p2c_call
+    {
+      LANG: Fortran
+    } "Testing conservative to primitive to conservative - MHD version"
+  } else {
+    schedule c2p2c in HydroBase_Initial AS c2p2c_call
+    {
+      LANG: Fortran
+    } "Testing conservative to primitive to conservative"
+  }
+}
+
+if (CCTK_Equals(initial_data,"prim2con2prim_test")) {
+  STORAGE:GRHydro_init_data_reflevel
+  schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE p2c2p_call
   {
     LANG: Fortran
-  } "Testing conservative to primitive to conservative"
+  } "Calculate current refinement level"
+
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
+  {
+    schedule p2c2pM in HydroBase_Initial AS p2c2p_call
+    {
+      LANG: Fortran
+    } "Testing primitive to conservative to primitive - MHD version"
+  } else {
+    schedule p2c2p in HydroBase_Initial AS p2c2p_call
+    {
+      LANG: Fortran
+    } "Testing primitive to conservative to primitive"
+  }
 }
 
 if (CCTK_Equals(initial_data,"reconstruction_test")) {
diff --git a/src/CheckParam.c b/src/CheckParam.c
index e758f4d..31222eb 100644
--- a/src/CheckParam.c
+++ b/src/CheckParam.c
@@ -11,5 +11,18 @@ void GRHydro_InitData_CheckParameters(CCTK_ARGUMENTS)
   {
       CCTK_PARAMWARN("You have to set 'HydroBase::timelevels to at least 2");
   }
+
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro") && 
+     ((CCTK_Equals(initial_hydro,"ony_atmo")) ||
+      (CCTK_Equals(initial_hydro,"read_conformal")) ||
+      (CCTK_Equals(initial_hydro,"simple_wave")) ||
+      (CCTK_Equals(initial_data,"con2primtest")) ||
+      (CCTK_Equals(initial_data,"reconstruction_test")) ||
+      (CCTK_Equals(shocktube_type,"diagshock")) ||
+      (CCTK_Equals(shocktube_type,"sphere")))) 
+    {
+      CCTK_PARAMWARN("That test not yet implemented in MHD!");
+    }
+	  
 }
 
diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90
new file mode 100644
index 0000000..12c1ed4
--- /dev/null
+++ b/src/GRHydro_C2P2CM.F90
@@ -0,0 +1,182 @@
+ /*@@
+   @file      GRHydro_C2P2CM.F90
+   @date      Sep 23, 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    c2p2cM
+   @date       Sep 23, 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 c2p2cM(CCTK_ARGUMENTS)
+
+  implicit none
+
+  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
+  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
+  CCTK_INT C2P_failed
+  logical epsnegative
+
+#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
+
+  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
+  
+  dens_send = 1.29047362d0
+  sx_send   = 0.166666658d0
+  sy_send   = 0.166666658d0
+  sz_send   = 0.166666658d0
+  tau_send  = 0.484123939d0
+  
+  bvcx_send   = Bx_init
+  bvcy_send   = By_init
+  bvcz_send   = Bz_init
+
+  eps_send = 1.0d-6
+  press_send = 6.666666666666667d-7
+  w_lorentz_send = 1.0d0
+
+  epsnegative = .false.
+  
+  GRHydro_rho_min = 1.e-10
+#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_EpsFromPress(GRHydro_eos_handle,keytemp,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(*,*) 'C2P2CM 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(*,*) 'C2P2CM test: getting the associated primitive variables.'
+  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,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,&
+       uxx,uxy,uxz,uyy,uyz,uzz,det,&
+       bvcx_send,bvcy_send,bvcz_send,b2_send,&
+       epsnegative,C2P_failed)
+  
+  write(*,*) 'C2P2CM 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(*,*) 'C2P2CM test: converting back to conserved variables.'
+  call prim2conM(GRHydro_eos_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(*,*) 'C2P2CM test: the conserved variables are'
+  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
+  
+  STOP
+
+  return
+
+end subroutine c2p2cM
diff --git a/src/GRHydro_Macros.h b/src/GRHydro_Macros.h
new file mode 100644
index 0000000..7009b86
--- /dev/null
+++ b/src/GRHydro_Macros.h
@@ -0,0 +1,12 @@
+#define SPATIAL_DETERMINANT(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_) \
+  (-(gxz_)**2*(gyy_) + 2*(gxy_)*(gxz_)*(gyz_) - (gxx_)*(gyz_)**2 - (gxy_)**2*(gzz_) \
+   + (gxx_)*(gyy_)*(gzz_))
+
+#define DOTP(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x1_,y1_,z1_,x2_,y2_,z2_) \
+ ( (gxx_)*(x1_)*(x2_)+(gyy_)*(y1_)*(y2_)+(gzz_)*(z1_)*(z2_)+ \
+   (gxy_)*( (x1_)*(y2_)+(y1_)*(x2_) )+(gxz_)*( (x1_)*(z2_)+(z1_)*(x2_) )+\
+   (gyz_)*( (y1_)*(z2_)+(z1_)*(y2_) ) )
+
+#define DOTP2(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x_,y_,z_)	\
+ ( (gxx_)*(x_)**2+(gyy_)*(y_)**2+(gzz_)*(z_)**2+ \
+  2.0*( (gxy_)*(x_)*(y_)+(gxz_)*(x_)*(z_)+(gyz_)*(y_)*(z_) ) )
diff --git a/src/GRHydro_P2C2P.F90 b/src/GRHydro_P2C2P.F90
new file mode 100644
index 0000000..10258c3
--- /dev/null
+++ b/src/GRHydro_P2C2P.F90
@@ -0,0 +1,167 @@
+ /*@@
+   @file      GRHydro_P2C2P.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    p2c2p
+   @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 p2c2p(CCTK_ARGUMENTS)
+
+  implicit none
+
+  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
+  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 pmin, epsmin
+  CCTK_INT C2P_failed
+  logical epsnegative
+
+#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
+
+  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
+  
+  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
+#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,&
+       GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,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'
+  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(*,*) 'P2C2P test: converting back to conserved variables.'
+  call prim2con(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+       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) 
+
+  write(*,*) 'P2C2P 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(*,*) 'P2C2P test: getting the associated primitive variables.'
+  call Con2Prim_pt(GRHydro_eos_handle,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, &
+       uxx,uxy,uxz,uyy,uyz,uzz,det,x_send,y_send,z_send,r_send,&
+       epsnegative,GRHydro_rho_min, pmin, epsmin, GRHydro_init_data_reflevel,C2P_failed)
+  
+  write(*,*) 'P2C2P 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
+  
+  STOP
+
+  return
+
+end subroutine p2c2p
diff --git a/src/GRHydro_P2C2PM.F90 b/src/GRHydro_P2C2PM.F90
new file mode 100644
index 0000000..c68b735
--- /dev/null
+++ b/src/GRHydro_P2C2PM.F90
@@ -0,0 +1,183 @@
+ /*@@
+   @file      GRHydro_P2C2PM.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(CCTK_ARGUMENTS)
+
+  implicit none
+
+  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
+  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
+  CCTK_INT C2P_failed
+  logical epsnegative
+
+#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
+
+  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
+#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,&
+       GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+  call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,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(*,*) '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 prim2conM(GRHydro_eos_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(*,*) 'P2C2PM test: getting the associated primitive variables.'
+  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,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,&
+       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
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90
new file mode 100644
index 0000000..1169fef
--- /dev/null
+++ b/src/GRHydro_ShockTubeM.F90
@@ -0,0 +1,174 @@
+ /*@@
+   @file      GRHydro_ShockTubeM.F90
+   @date      Sep 23, 2010
+   @author    Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+   @desc 
+   Initial data of the shock tube type - MHD version.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#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)
+#define velz(i,j,k) vel(i,j,k,3)
+#define sx(i,j,k) scon(i,j,k,1)
+#define sy(i,j,k) scon(i,j,k,2)
+#define sz(i,j,k) scon(i,j,k,3)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
+
+ /*@@
+   @routine    GRHydro_shocktubeM
+   @date       Sat Jan 26 02:53:49 2002
+   @author     Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+   @desc 
+   Initial data for shock tubes, parallel to
+   a coordinate axis. Either Sods problem or the standard shock tube.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Expansion and alteration of the test code from GRAstro_Hydro, 
+   written by Mark Miller.
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  CCTK_INT :: i, j, k, nx, ny, nz
+  CCTK_REAL :: direction, det
+  CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, &
+       velzl, velzr, epsl, epsr
+  CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  do i=1,nx
+    do j=1,ny
+      do k=1,nz
+        if (CCTK_EQUALS(shocktube_type,"diagshock")) then
+          direction = x(i,j,k) - shock_xpos + &
+               y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos
+        else if (CCTK_EQUALS(shocktube_type,"xshock")) then
+          direction = x(i,j,k) - shock_xpos
+        else if (CCTK_EQUALS(shocktube_type,"yshock")) then
+          direction = y(i,j,k) - shock_ypos
+        else if (CCTK_EQUALS(shocktube_type,"zshock")) then
+          direction = z(i,j,k) - shock_zpos
+        else if (CCTK_EQUALS(shocktube_type,"sphere")) then
+          direction = sqrt((x(i,j,k)-shock_xpos)**2+&
+                           (y(i,j,k)-shock_ypos)**2+&
+                           (z(i,j,k)-shock_zpos)**2)-shock_radius
+        end if
+        if (CCTK_EQUALS(shock_case,"Simple")) then
+          rhol = 10.d0
+          rhor = 1.d0
+          velxl = 0.d0
+          velxr = 0.d0
+          velyl = 0.d0
+          velyr = 0.d0
+          velzl = 0.d0
+          velzr = 0.d0
+          epsl = 2.d0
+          epsr = 1.d-6
+        else if (CCTK_EQUALS(shock_case,"Sod")) 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
+          epsl = 1.5d0
+          epsr = 1.2d0
+!!$This line only for polytrope, k=1
+!!$          epsr = 0.375d0
+        else if (CCTK_EQUALS(shock_case,"Blast")) 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
+          epsl = 1500.d0
+          epsr = 1.5d-2
+        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
+          velx(i,j,k) = velxl
+          vely(i,j,k) = velyl
+          velz(i,j,k) = velzl
+          eps(i,j,k) = epsl
+          Bvecx(i,j,k)=bvcxl
+          Bvecy(i,j,k)=bvcyl
+          Bvecz(i,j,k)=bvczl
+        else
+          rho(i,j,k) = rhor
+          velx(i,j,k) = velxr
+          vely(i,j,k) = velyr
+          velz(i,j,k) = velzr
+          eps(i,j,k) = epsr
+          Bvecx(i,j,k)=bvcxr
+          Bvecy(i,j,k)=bvcyr
+          Bvecz(i,j,k)=bvczr
+        end if
+
+        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 Prim2ConPolyM(GRHydro_eos_handle,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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+               tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),&
+               velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+               eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+        else
+          call Prim2ConGenM(GRHydro_eos_handle,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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+               tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),&
+               velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+               eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+        end if
+    enddo
+    enddo
+  enddo
+
+  densrhs = 0.d0
+  srhs = 0.d0
+  taurhs = 0.d0
+  Bvecrhs = 0.d0
+
+
+  return
+  
+end subroutine GRHydro_shocktubeM
diff --git a/src/make.code.defn b/src/make.code.defn
index fd18296..dc78381 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -3,13 +3,17 @@
 
 # Source files in this directory
 SRCS = 	GRHydro_C2P2C.F90 \
+	GRHydro_C2P2CM.F90 \
 	GRHydro_Con2Prim.F90 \
 	GRHydro_ReconstructTest.F90 \
 	GRHydro_ShockTube.F90 \
+	GRHydro_ShockTubeM.F90 \
 	GRHydro_Only_Atmo.F90 \
 	GRHydro_ReadConformalData.F90 \
 	GRHydro_SimpleWave.F90 \
-  CheckParam.c
+   CheckParam.c\
+   GRHydro_P2C2P.F90 \
+   GRHydro_P2C2PM.F90
 
 # Subdirectories containing source files
 SUBDIRS =