RIT MHD dev:

 Balsara's and Komissarov's shock tube tests.


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

7 files changed, 657 insertions, 23 deletions

diff --git a/par/balsara1a.par b/par/balsara1a.par
index 1943166..6997407 100644
--- a/par/balsara1a.par
+++ b/par/balsara1a.par
@@ -5,7 +5,8 @@ ActiveThorns = "time MoL cartgrid3d carpetioascii ioutil Fortran boundary  hydro
 driver::ghost_size=3
 grhydro::grhydro_stencil=3
 
-time::dtfac = 0.25
+#time::dtfac = 0.25
+time::dtfac = 0.8
 
 methodoflines::ODE_Method = "rk2"
 methodoflines::MoL_Intermediate_Steps=2
@@ -37,25 +38,30 @@ grhydro::GRHydro_MaxNumEvolvedVars=10
 
 grid::type = "BySpacing"
 grid::domain = "full"
-grid::dxyz = 0.01
+#grid::dxyz = 0.01
+grid::dxyz = 0.000625
 
-driver::global_nx = 150
+driver::global_nx = 1600
 driver::global_ny = 10
 driver::global_nz = 10
 
+Cactus::terminate="time"
+Cactus::cctk_final_time = 0.4
 #cactus::cctk_itlast = 40
-cactus::cctk_itlast = 200
+#cactus::cctk_itlast = 200
+#cactus::cctk_itlast = 2
 
-IO::out_dir = "balsara1a"
+IO::out_dir = $parfile
 #IOBasic::outInfo_every = 1
 #IOBasic::outInfo_vars = "HydroBase::rho"
 CarpetIOBasic::outInfo_vars="hydrobase::rho"
 CarpetIOBasic::outInfo_every=1
 CarpetIOASCII::out1D_criterion = "time"
 CarpetIOASCII::out1D_dt = 0.01
+CarpetIOASCII::out1D_d=no
 CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec"
 
 CarpetIOHDF5::out_every = 1
-CarpetIOHDF5::out_vars  = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec"
+CarpetIOHDF5::out_vars  = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel HydroBase::w_lorentz grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec"
 
 
diff --git a/param.ccl b/param.ccl
index 01fbe7f..8d4da47 100644
--- a/param.ccl
+++ b/param.ccl
@@ -12,6 +12,8 @@ EXTENDS KEYWORD initial_hydro ""
   "only_atmo"	:: "Set only a low atmosphere"
   "read_conformal":: "After reading in initial alp, rho and gxx from h5 files, sets the other quantities"
   "simple_wave" :: "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)"
+  "Monopole"    :: "Monopole at the center"
+  "cylexp"      :: "Cylindrical Explosion"
 }
 
 shares:ADMBase
@@ -31,22 +33,35 @@ private:
 KEYWORD shocktube_type "Diagonal or parallel shock?"
 {
   "diagshock"	:: "Diagonal across all axes"
+  "diagshock2d"	:: "Diagonal across x-y axes"
   "xshock"	:: "Parallel to x axis"
   "yshock"	:: "Parallel to y axis"
   "zshock"	:: "Parallel to z axis"
   "sphere"  :: "spherically symmetric shock"
-} "diagshock"
+} "xshock"
 
 KEYWORD shock_case "Simple, Sod's problem or other?"
 {
   "Simple"	:: "GRAstro_Hydro test case"
   "Sod"		:: "Sod's problem"
   "Blast"	:: "Strong blast wave"
+  "Balsaralike1" :: "Hydro version of Balsara Test #1"
+  "Balsara0"    :: "Balsara Test #1, but unmagnetized"
   "Balsara1"    :: "Balsara Test #1"
   "Balsara2"    :: "Balsara Test #2"
   "Balsara3"    :: "Balsara Test #3"
   "Balsara4"    :: "Balsara Test #4"
   "Balsara5"    :: "Balsara Test #5"
+  "Alfven"      :: "Generical Alfven Test"
+  "Komissarov1" :: "Komissarov Test #1"
+  "Komissarov2" :: "Komissarov Test #2"
+  "Komissarov3" :: "Komissarov Test #3"
+  "Komissarov4" :: "Komissarov Test #4"
+  "Komissarov5" :: "Komissarov Test #5"
+  "Komissarov6" :: "Komissarov Test #6"
+  "Komissarov7" :: "Komissarov Test #7"
+  "Komissarov8" :: "Komissarov Test #8"
+  "Komissarov9" :: "Komissarov Test #9"
 } "Sod"
 
 REAL shock_xpos "Position of shock plane: x"
@@ -124,3 +139,5 @@ USES real rho_rel_min
 USES REAL initial_rho_abs_min
 USES REAL initial_rho_rel_min
 USES REAL initial_atmosphere_factor
+USES int GRHydro_stencil
+USES BOOLEAN clean_divergence
diff --git a/schedule.ccl b/schedule.ccl
index a35dfdb..7040982 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -5,6 +5,13 @@ schedule GRHydro_InitData_CheckParameters  AT CCTK_PARAMCHECK
   LANG: C
 } "Check parameters"
 
+if (CCTK_Equals(initial_hydro,"Monopole")) {
+  schedule GRHydro_MonopoleM in HydroBase_Initial
+    {
+      LANG: Fortran
+    } "Monopole initial data - MHD version"
+}
+
 if (CCTK_Equals(initial_hydro,"shocktube")) {
   if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
   {
@@ -13,6 +20,22 @@ if (CCTK_Equals(initial_hydro,"shocktube")) {
       LANG: Fortran
     } "Shocktube initial data - MHD version"
 
+    if(CCTK_Equals(shocktube_type,"diagshock")) 
+    {
+      schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+      {
+        LANG: Fortran	
+      } "Diagonal shock boundary conditions"
+    }
+
+    if(CCTK_Equals(shocktube_type,"diagshock2d")) 
+    {
+      schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+      {
+        LANG: Fortran	
+      } "2-D Diagonal shock boundary conditions"
+    }
+
   } else {
     schedule GRHydro_shocktube in HydroBase_Initial
     {
@@ -22,6 +45,13 @@ if (CCTK_Equals(initial_hydro,"shocktube")) {
   }
 }
 
+if (CCTK_Equals(initial_hydro,"cylexp")) {
+  schedule GRHydro_CylindricalExplosionM in HydroBase_Initial
+    {
+      LANG: Fortran
+    } "Cylindrical Explosion initial data - MHD-only"
+}
+
 if (CCTK_Equals(initial_data,"con2primtest")) {
   STORAGE:GRHydro_init_data_reflevel
   schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE GRHydro_con2primtest
diff --git a/src/CheckParam.c b/src/CheckParam.c
index 124db7f..7182dab 100644
--- a/src/CheckParam.c
+++ b/src/CheckParam.c
@@ -18,20 +18,31 @@ void GRHydro_InitData_CheckParameters(CCTK_ARGUMENTS)
       (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!");
     }
 
   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_Equals(shock_case,"Balsara0"   )) ||
+      (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_Equals(shock_case,"Alfven"     )) ||
+      (CCTK_Equals(shock_case,"Komissarov1")) ||
+      (CCTK_Equals(shock_case,"Komissarov2")) ||
+      (CCTK_Equals(shock_case,"Komissarov3")) ||
+      (CCTK_Equals(shock_case,"Komissarov4")) ||
+      (CCTK_Equals(shock_case,"Komissarov5")) ||
+      (CCTK_Equals(shock_case,"Komissarov6")) ||
+      (CCTK_Equals(shock_case,"Komissarov7")) ||
+      (CCTK_Equals(shock_case,"Komissarov8")) ||
+      (CCTK_Equals(shock_case,"Komissarov9")) ||
+      (CCTK_Equals(initial_hydro,"cylexp")) 
+      ))
     {
       CCTK_PARAMWARN("That test requires MHD!  Set Bvec_evolution_method=GRHYDRO!");
     }
 }
-
diff --git a/src/GRHydro_ShockTube.F90 b/src/GRHydro_ShockTube.F90
index baf891d..8b9105b 100644
--- a/src/GRHydro_ShockTube.F90
+++ b/src/GRHydro_ShockTube.F90
@@ -93,6 +93,17 @@ subroutine GRHydro_shocktube(CCTK_ARGUMENTS)
           velzr = 0.d0
           epsl = 1.5d0
           epsr = 1.2d0
+        else if (CCTK_EQUALS(shock_case,"Balsaralike1")) then
+          rhol = 1.0d0
+          rhor = 0.125d0
+          velxl = 0.0d0
+          velxr = 0.0d0
+          velyl = 0.0d0
+          velyr = 0.0d0
+          velzl = 0.0d0
+          velzr = 0.0d0
+          epsl = 1.0d0/rhol
+          epsr = 0.1d0/rhor
 !!$This line only for polytrope, k=1
 !!$          epsr = 0.375d0
         else if (CCTK_EQUALS(shock_case,"Blast")) then
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90
index 186cf69..72bac3e 100644
--- a/src/GRHydro_ShockTubeM.F90
+++ b/src/GRHydro_ShockTubeM.F90
@@ -23,6 +23,11 @@
 #define Bvecy(i,j,k) Bvec(i,j,k,2)
 #define Bvecz(i,j,k) Bvec(i,j,k,3)
 
+#define OOSQRT2 (0.7071067811865475244008442)
+#define OOSQRT3 (0.5773502691896257645091489)
+#define OOSQRT6 (0.4082482904638630163662140)
+
+
  /*@@
    @routine    GRHydro_shocktubeM
    @date       Sat Jan 26 02:53:49 2002
@@ -53,6 +58,8 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
   CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, &
        velzl, velzr, epsl, epsr
   CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr
+  CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3
+  
   
   bvcxl = Bx_init
   bvcyl = By_init
@@ -69,18 +76,27 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
     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
+!!$ The diagshock choice yields a shock plane perpendicular to the fixed vector (1,1,1)
+!!$ This could be changed, but would require 3 new params containing the new shock direction
+           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,"diagshock2d")) then
+!!$ The diagshock choice yields a shock plane perpendicular to the fixed vector (1,1,0), with similarity in the z-dir.
+!!$ This could be changed, but would require 2 new params containing the new shock direction
+           direction = x(i,j,k) - shock_xpos + &
+                y(i,j,k) - shock_ypos
+
         else if (CCTK_EQUALS(shocktube_type,"xshock")) then
-          direction = x(i,j,k) - shock_xpos
+           direction = x(i,j,k) - shock_xpos
         else if (CCTK_EQUALS(shocktube_type,"yshock")) then
-          direction = y(i,j,k) - shock_ypos
+           direction = y(i,j,k) - shock_ypos
         else if (CCTK_EQUALS(shocktube_type,"zshock")) then
-          direction = z(i,j,k) - shock_zpos
+           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
+           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
@@ -117,6 +133,31 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           velzr = 0.d0
           epsl = 1500.d0
           epsr = 1.5d-2
+
+!!$ The following shocktubes are from Balsara 2001 .
+!!$   All use n=1600 cells, over domain x=[-0.5,0.5]
+!!$   All assume ideal-gas or gamma-law EOS, the first test uses GAMMA=2. while the rest GAMMA=5./3.
+
+!!$ Unmagnetized Test 1 (rel. Brio & Wu 1988 by van Putten 1993) of Balsara 2001 -- compare at t=0.4
+       else if (CCTK_EQUALS(shock_case,"Balsara0")) then
+          rhol = 1.0d0
+          rhor = 0.125d0
+          velxl = 0.0d0
+          velxr = 0.0d0
+          velyl = 0.0d0
+          velyr = 0.0d0
+          velzl = 0.0d0
+          velzr = 0.0d0
+          bvcxl=0.0d0
+          bvcxr=0.0d0
+          bvcyl=0.0d0
+          bvcyr=0.0d0
+          bvczl=0.0d0
+          bvczr=0.0d0
+          epsl = 1.0d0/rhol
+          epsr = 0.1d0/rhor
+
+!!$ Test 1 (rel. Brio & Wu 1988 by van Putten 1993) of Balsara 2001 -- compare at t=0.4
        else if (CCTK_EQUALS(shock_case,"Balsara1")) then
           rhol = 1.0d0
           rhor = 0.125d0
@@ -134,6 +175,8 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           bvczr=0.0d0
           epsl = 1.0d0/rhol
           epsr = 0.1d0/rhor
+
+!!$ Test 2 (blast wave) of Balsara 2001  -- compare at t=0.4
        else if (CCTK_EQUALS(shock_case,"Balsara2")) then
           rhol = 1.d0
           rhor = 1.d0
@@ -151,6 +194,8 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           bvczr=0.7d0
           epsl = 1.5d0*30.0d0/rhol
           epsr = 1.5d0*1.0d0/rhor
+
+!!$ Test 3 (blast wave) of Balsara 2001  -- compare at t=0.4
        else if (CCTK_EQUALS(shock_case,"Balsara3")) then
           rhol = 1.d0
           rhor = 1.d0
@@ -168,6 +213,8 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           bvczr=0.7d0
           epsl = 1.5d0*1000.0d0/rhol
           epsr = 1.5d0*0.1d0/rhor
+
+!!$ Test 4 (rel. version of Noh 1987) of Balsara 2001  -- compare at t=0.4
        else if (CCTK_EQUALS(shock_case,"Balsara4")) then
           rhol = 1.d0
           rhor = 1.d0
@@ -185,6 +232,8 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           bvczr=-7.d0
           epsl = 1.5d0*0.1d0/rhol
           epsr = 1.5d0*0.1d0/rhor
+
+!!$ Test 5 (non-coplanar set of waves) of Balsara 2001  -- compare at t=0.55
        else if (CCTK_EQUALS(shock_case,"Balsara5")) then
           rhol = 1.08d0
           rhor = 1.d0
@@ -202,12 +251,286 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           bvczr=0.5d0
           epsl = 1.5d0*0.95d0/rhol
           epsr = 1.5d0*1.0d0/rhor
+
+!!$  "Generic Alfven Test of  Giacomazzo and Rezzolla J.Comp.Phys (2006)
+       else if (CCTK_EQUALS(shock_case,"Alfven")) then
+          rhol = 1.d0
+          rhor = 0.9d0
+          velxl = 0.d0
+          velxr = 0.d0
+          velyl = 0.3d0
+          velyr = 0.d0
+          velzl = 0.4d0
+          velzr = 0.d0
+          bvcxl=1.0d0
+          bvcxr=1.0d0
+          bvcyl=6.d0
+          bvcyr=5.d0
+          bvczl=2.d0
+          bvczr=2.d0
+          epsl = 1.5d0*5.d0/rhol
+          epsr = 1.5d0*5.3d0/rhor
+
+!!$ The following 9 tests are from Komissarov 1999 . 
+!!$   Note that the data is specified in terms of the 4-velocity, so some conversion is necessary. 
+!!$   All assume ideal-gas or gamma-law EOS with GAMMA=4./3. 
+
+!!$ Fast Shock test of Komissarov 1999 -- compare  at t=2.5 , n=40, x=[-1,1]
+       else if (CCTK_EQUALS(shock_case,"Komissarov1")) then
+          rhol = 1.d0
+          rhor = 25.48d0
+          bvcxl=20.d0
+          bvcxr=20.d0
+          bvcyl=25.02d0
+          bvcyr=49.d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d0  /rhol
+          epsr = 3.d0 * 367.5d0 /rhor
+
+          ux=25.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=1.091d0
+          uy=0.3923d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
+!!$ Slow Shock test of Komissarov 1999 -- compare at t=2. , n=200,  x=[-0.5,1.5]
+       else if (CCTK_EQUALS(shock_case,"Komissarov2")) then
+          rhol = 1.d0
+          rhor = 3.323d0
+          bvcxl=10.d0
+          bvcxr=10.d0
+          bvcyl=18.28d0
+          bvcyr=14.49d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d0  /rhol
+          epsr = 3.d0 * 55.36d0 /rhor
+
+          ux=1.53d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=0.9571d0
+          uy=-0.6822d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
+!!$ Switch-off Fast test of Komissarov 1999 -- compare  at t=1. , n=150, x=[-1,1] 
+       else if (CCTK_EQUALS(shock_case,"Komissarov3")) then
+          rhol = 0.1d0
+          rhor = 0.562d0
+          bvcxl=2.d0
+          bvcxr=2.d0
+          bvcyl=0.d0
+          bvcyr=4.710d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d0  /rhol
+          epsr = 3.d0 * 10.d0 /rhor
+
+          ux=-2.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=-0.212d0
+          uy=-0.590d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
+!!$ Switch-on Slow test of Komissarov 1999 -- compare  at t=2. , n=150, x=[-1,1.5]
+       else if (CCTK_EQUALS(shock_case,"Komissarov4")) then
+          rhol = 1.78d-2
+          rhor = 1.d-2
+          bvcxl=1.d0
+          bvcxr=1.d0
+          bvcyl=1.022d0
+          bvcyr=0.d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 0.1d0  /rhol
+          epsr = 3.d0 * 1.d0 /rhor
+
+          ux=-0.765d0
+          uy=-1.386d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=0.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
+!!$ Alfven wave test of Komissarov 1999 -- compare  at t=2. , n=200, x=[-1,1.5]
+!!$   Needs special setup -- FIX
+       else if (CCTK_EQUALS(shock_case,"Komissarov5")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          bvcxl=3.d0
+          bvcxr=3.d0
+          bvcyl=3.d0
+          bvcyr=-6.857d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d0  /rhol
+          epsr = 3.d0 * 1.d0 /rhor
+
+          ux=0.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=3.70d0
+          uy=5.76d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
+!!$ Compound wave test of Komissarov 1999 -- compare  at t=0.1,0.75,1.5 , n=200, x=[-0.5,1.5]
+!!$   Needs special setup -- FIX
+       else if (CCTK_EQUALS(shock_case,"Komissarov6")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          bvcxl=3.d0
+          bvcxr=3.d0
+          bvcyl=3.d0
+          bvcyr=-6.857d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d0  /rhol
+          epsr = 3.d0 * 1.d0 /rhor
+
+          ux=0.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=3.70d0
+          uy=5.76d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
+!!$ Shock Tube 1  test of Komissarov 1999 -- compare  at t=1. , n=400, x=[-1,1.5]
+       else if (CCTK_EQUALS(shock_case,"Komissarov7")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          bvcxl=1.d0
+          bvcxr=1.d0
+          bvcyl=0.d0
+          bvcyr=0.d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d3  /rhol
+          epsr = 3.d0 * 1.d0 /rhor
+
+          velxl = 0.
+          velyl = 0.
+          velzl = 0.
+
+          velxr = 0.
+          velyr = 0.
+          velzr = 0.
+
+!!$ Shock Tube 2 test of Komissarov 1999 -- compare  at t=1. , n=500, x=[-1.25,1.25]
+       else if (CCTK_EQUALS(shock_case,"Komissarov8")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          bvcxl=0.d0
+          bvcxr=0.d0
+          bvcyl=20.d0
+          bvcyr=0.d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 30.d0  /rhol
+          epsr = 3.d0 * 1.d0 /rhor
+
+          velxl = 0.
+          velyl = 0.
+          velzl = 0.
+
+          velxr = 0.
+          velyr = 0.
+          velzr = 0.
+
+!!$ Collision test of Komissarov 1999 -- compare  at t=1.22 , n=200, x=[-1,1]
+       else if (CCTK_EQUALS(shock_case,"Komissarov9")) then
+          rhol = 1.d0
+          rhor = 1.d0
+          bvcxl=10.d0
+          bvcxr=10.d0
+          bvcyl=10.d0
+          bvcyr=-10.d0
+          bvczl=0.d0
+          bvczr=0.d0
+          epsl = 3.d0 * 1.d0  /rhol
+          epsr = 3.d0 * 1.d0 /rhor
+
+          ux=5.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxl = ux/ut
+          velyl = uy/ut
+          velzl = uz/ut
+
+          ux=-5.d0
+          uy=0.d0
+          uz=0.d0
+          ut = sqrt(1. + ux*ux + uy*uy + uz*uz)
+          velxr = ux/ut
+          velyr = uy/ut
+          velzr = uz/ut
+
        else
           call CCTK_WARN(0,"Shock case not recognized")
         end if
 
         if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& 
              ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then
+
+!!$ Left state
+
           rho(i,j,k) = rhol
           velx(i,j,k) = velxl
           vely(i,j,k) = velyl
@@ -217,6 +540,9 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           Bvecy(i,j,k)=bvcyl
           Bvecz(i,j,k)=bvczl
         else
+
+!!$ Right state
+
           rho(i,j,k) = rhor
           velx(i,j,k) = velxr
           vely(i,j,k) = velyr
@@ -227,6 +553,58 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           Bvecz(i,j,k)=bvczr
         end if
 
+        
+        if (CCTK_EQUALS(shocktube_type,"yshock")) then
+!!$ Cycle x,y,z forward           
+           tmp=velx(i,j,k)
+           velx(i,j,k)=velz(i,j,k)
+           velz(i,j,k)=vely(i,j,k)
+           vely(i,j,k)=tmp
+           tmp=Bvecx(i,j,k)
+           Bvecx(i,j,k)=Bvecz(i,j,k)
+           Bvecz(i,j,k)=Bvecy(i,j,k)
+           Bvecy(i,j,k)=tmp
+        else if (CCTK_EQUALS(shocktube_type,"zshock")) then
+!!$ Cycle x,y,z backward           
+           tmp=velx(i,j,k)
+           velx(i,j,k)=vely(i,j,k)
+           vely(i,j,k)=velz(i,j,k)
+           velz(i,j,k)=tmp
+           tmp=Bvecx(i,j,k)
+           Bvecx(i,j,k)=Bvecy(i,j,k)
+           Bvecy(i,j,k)=Bvecz(i,j,k)
+           Bvecz(i,j,k)=tmp
+        else if (CCTK_EQUALS(shocktube_type,"diagshock")) then
+!!$ Rotated basis vectors necessary to evaluate the orthogonal matrix elements:
+!!$ xhat = 1/sqrt(3)[1,1,1], yhat = 1/sqrt(2)[-1,1,0]; zhat = 1/sqrt(6)[-1,-1,2]
+!!$ Orthogonal matrix constructed from the tensor product between the original
+!!$ cartesian basis x's and new basis vectors xhat's, rotated towards the diagonal 
+!!$ shock normal and tangent directions.
+           tmp  = OOSQRT3*velx(i,j,k) - OOSQRT2*vely(i,j,k) - OOSQRT6*velz(i,j,k)
+           tmp2 = OOSQRT3*velx(i,j,k) + OOSQRT2*vely(i,j,k) - OOSQRT6*velz(i,j,k)
+           tmp3 = OOSQRT3*velx(i,j,k)                  + 2.d0*OOSQRT6*velz(i,j,k)
+           velx(i,j,k)=tmp
+           vely(i,j,k)=tmp2
+           velz(i,j,k)=tmp3
+           tmp  = OOSQRT3*Bvecx(i,j,k) - OOSQRT2*Bvecy(i,j,k) - OOSQRT6*Bvecz(i,j,k)
+           tmp2 = OOSQRT3*Bvecx(i,j,k) + OOSQRT2*Bvecy(i,j,k) - OOSQRT6*Bvecz(i,j,k)
+           tmp3 = OOSQRT3*Bvecx(i,j,k)                   + 2.d0*OOSQRT6*Bvecz(i,j,k)
+           Bvecx(i,j,k)=tmp
+           Bvecy(i,j,k)=tmp2
+           Bvecz(i,j,k)=tmp3
+        else if (CCTK_EQUALS(shocktube_type,"diagshock2d")) then
+!!$ New basis:
+!!$ xhat = 1/sqrt(2)[1,1,0], yhat = 1/sqrt(2)[-1,1,0]; zhat = [0,0,1]
+           tmp  = OOSQRT2*velx(i,j,k) - OOSQRT2*vely(i,j,k)
+           tmp2 = OOSQRT2*velx(i,j,k) + OOSQRT2*vely(i,j,k)
+           velx(i,j,k)=tmp
+           vely(i,j,k)=tmp2
+           tmp  = OOSQRT2*Bvecx(i,j,k) - OOSQRT2*Bvecy(i,j,k)
+           tmp2 = OOSQRT2*Bvecx(i,j,k) + OOSQRT2*Bvecy(i,j,k)
+           Bvecx(i,j,k)=tmp
+           Bvecy(i,j,k)=tmp2
+        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))
 
         if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then
@@ -257,3 +635,182 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
   return
   
 end subroutine GRHydro_shocktubeM
+
+subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  CCTK_INT :: i, j, k, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew
+  CCTK_INT :: xoff,yoff,zoff,indsum
+  CCTK_REAL :: det
+
+  stenp1=GRHydro_stencil + 1
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  minsum = 3*stenp1
+  maxsum = nx+ny+nz-3*stenp1 + 3
+
+  xoff=0
+  yoff=0
+  zoff=0
+        
+  do k=1,nz
+     if(k.lt.stenp1)then
+              zoff=k-stenp1
+           else if(k.gt.nz-stenp1+1) then
+              zoff=k-(nz-stenp1+1)
+           else
+              zoff=0
+           endif
+
+     do j=1,ny
+        if(j.lt.stenp1)then
+           yoff=j-stenp1
+        else if(j.gt.ny-stenp1+1) then
+           yoff=j-(ny-stenp1+1)
+        else
+           yoff=0
+        endif
+
+        do i=1,nx
+           if(i.lt.stenp1) then
+              xoff=i-stenp1
+           else if(i.gt.nx-stenp1+1) then
+              xoff=i-(nx-stenp1+1)
+           else
+              xoff=0
+           endif
+
+           indsum = i+j+k
+           
+           if( (xoff.ne.0.or.yoff.ne.0.or.zoff.ne.0) .and. &
+                indsum.ge.minsum.and.indsum.le.maxsum) then
+!!$ We can map the point to the interior diagonal, orthogonal to the shock.
+
+              inew=indsum/3
+              jnew=(indsum-inew)/2
+              knew=indsum-inew-jnew
+              
+              dens(i,j,k) = dens(inew,jnew,knew)
+              sx(i,j,k) = sx(inew,jnew,knew)
+              sy(i,j,k) = sy(inew,jnew,knew)
+              sz(i,j,k) = sz(inew,jnew,knew)
+              tau(i,j,k) = tau(inew,jnew,knew)
+              Bvecx(i,j,k)=Bvecx(inew,jnew,knew)
+              Bvecy(i,j,k)=Bvecy(inew,jnew,knew)
+              Bvecz(i,j,k)=Bvecz(inew,jnew,knew)
+             if(clean_divergence.ne.0) then 
+               psidc(i,j,k)=psidc(inew,jnew,knew)
+             endif
+
+           endif
+           
+        enddo
+     enddo
+  enddo
+  
+end subroutine GRHydro_Diagshock_BoundaryM
+
+subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  CCTK_INT :: i, j, k, kc, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew
+  CCTK_INT :: xoff,yoff,zoff,indsum
+  CCTK_REAL :: det
+  
+  stenp1=GRHydro_stencil + 1
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  minsum = 2*stenp1
+  maxsum = nx+ny-2*stenp1 + 2
+  
+  xoff=0
+  yoff=0
+  zoff=0
+        
+  do k=1,nz
+     if(k.lt.stenp1)then
+              zoff=k-stenp1
+           else if(k.gt.nz-stenp1+1) then
+              zoff=k-(nz-stenp1+1)
+           else
+              zoff=0
+           endif
+
+     do j=1,ny
+        if(j.lt.stenp1)then
+           yoff=j-stenp1
+        else if(j.gt.ny-stenp1+1) then
+           yoff=j-(ny-stenp1+1)
+        else
+           yoff=0
+        endif
+
+        do i=1,nx
+           if(i.lt.stenp1) then
+              xoff=i-stenp1
+           else if(i.gt.nx-stenp1+1) then
+              xoff=i-(nx-stenp1+1)
+           else
+              xoff=0
+           endif
+  
+        indsum = i+j
+        
+        if( (xoff.ne.0.or.yoff.ne.0) .and. &
+             indsum.ge.minsum.and.indsum.le.maxsum) then
+!!$ We can map the point to the interior diagonal, orthogonal to the shock.
+           
+           inew=indsum/2
+           jnew=indsum-inew
+           
+           dens(i,j,k) = dens(inew,jnew,k)
+           sx(i,j,k) = sx(inew,jnew,k)
+           sy(i,j,k) = sy(inew,jnew,k)
+           sz(i,j,k) = sz(inew,jnew,k)
+           tau(i,j,k) = tau(inew,jnew,k)
+           Bvecx(i,j,k)=Bvecx(inew,jnew,k)
+           Bvecy(i,j,k)=Bvecy(inew,jnew,k)
+           Bvecz(i,j,k)=Bvecz(inew,jnew,k)
+           if(clean_divergence.ne.0) then 
+              psidc(i,j,k)=psidc(inew,jnew,k)
+           endif
+
+        else if( zoff.ne.0) then
+ 
+           kc = nz/2
+
+           dens(i,j,k) = dens(i,j,kc)
+           sx(i,j,k) = sx(i,j,kc)
+           sy(i,j,k) = sy(i,j,kc)
+           sz(i,j,k) = sz(i,j,kc)
+           tau(i,j,k) = tau(i,j,kc)
+           Bvecx(i,j,k)=Bvecx(i,j,kc)
+           Bvecy(i,j,k)=Bvecy(i,j,kc)
+           Bvecz(i,j,k)=Bvecz(i,j,kc)
+           if(clean_divergence.ne.0) then
+              psidc(i,j,k)=psidc(i,j,kc)
+           endif
+           
+        endif
+        
+        enddo
+     enddo
+  enddo
+  
+end subroutine GRHydro_Diagshock2D_BoundaryM
diff --git a/src/make.code.defn b/src/make.code.defn
index 09a09f2..7bbf06c 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -14,7 +14,9 @@ SRCS = 	GRHydro_C2P2C.F90 \
    CheckParam.c\
    GRHydro_P2C2P.F90 \
    GRHydro_P2C2PM.F90 \
-   GRHydro_P2C2PM_polytype.F90
+   GRHydro_P2C2PM_polytype.F90 \
+   GRHydro_MonopoleM.F90 \
+   GRHydro_CylindricalExplosionM.F90
 
 # Subdirectories containing source files
 SUBDIRS =