1 files changed, 109 insertions, 0 deletions
diff --git a/src/GRHydro_MonopoleM.F90 b/src/GRHydro_MonopoleM.F90
new file mode 100644
index 0000000..a34a51c
--- /dev/null
+++ b/src/GRHydro_MonopoleM.F90
@@ -0,0 +1,109 @@
+ /*@@
+   @file      GRHydro_MonopoleM.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_MonopoleM
+   @date       Sat Jan 26 02:53:49 2002
+   @author     Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+   @desc 
+   A monopole in space
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Expansion and alteration of the test code from GRAstro_Hydro, 
+   written by Mark Miller.
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_MonopoleM(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
+  CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3
+  
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  do i=1,nx
+     do j=1,ny
+        do k=1,nz
+           
+           rho(i,j,k) = 1.0
+           velx(i,j,k) = 0.0
+           vely(i,j,k) = 0.0
+           velz(i,j,k) = 0.0
+           eps(i,j,k) = 0.1
+           if(i.eq.nx/2+1.and.j.eq.ny/2+1.and.k.eq.nz/2+1) then
+              Bvecx(i,j,k)=0.1
+           else
+              Bvecx(i,j,k)=0.0
+           endif
+           Bvecy(i,j,k)=0.0
+           Bvecz(i,j,k)=0.0
+
+           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_MonopoleM