RIT MHD dev:

1) Komissarov's cylindrical explosion test
2) Magnetic monopole test



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

2 files changed, 273 insertions, 0 deletions

diff --git a/src/GRHydro_CylindricalExplosionM.F90 b/src/GRHydro_CylindricalExplosionM.F90
new file mode 100644
index 0000000..d06f777
--- /dev/null
+++ b/src/GRHydro_CylindricalExplosionM.F90
@@ -0,0 +1,164 @@
+ /*@@
+   @file      GRHydro_CylindricalExplosionM.F90
+   @date      Apr 22, 2011
+   @author    Scott Noble, Joshua Faber, Bruno Mundim
+   @desc 
+   Cylindrical magnetized shocks.
+   @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_cylindricalexplosionM
+   @date       Fri Apr 22 14:51:37 EDT 2011
+   @author     Scott Noble, Joshua Faber, Bruno Mundim
+   @desc 
+   Initial data for cylidrically symmetric shocks with 
+   magnetic fields.  Meant to recreate tests demonstrated 
+   in Komissarov (1999).
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Using GRHydro_ShockTubeM.F90 as a template.
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_cylindricalexplosionM(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, pressl, pressr
+  CCTK_REAL :: bvcxl,bvcyl,bvczl
+  CCTK_REAL :: tmp,tmp2,gam,r_inner,r_outer
+
+  CCTK_REAL :: cyl_fr
+  
+!!$Uses Bx_init, By_init and Bz_init to set magnetic field strength
+!!$Original tests had Bx_init = 0.1, 1.0   with By_init=Bz_init = 0
+  
+  bvcxl = Bx_init
+  bvcyl = By_init
+  bvczl = Bz_init
+  
+!!$Adiabatic index for test
+  gam = (4.d0/3.d0)
+
+!!$Inner radius and outer radius
+  r_inner = 8.d-1
+  r_outer = 1.d0
+  
+!!$Inner values
+  rhol   = 1.d-4
+  pressl = 3.d-5
+
+!!$Outer values
+  rhor   = 1.d-2
+  pressr = 1.d0
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+
+  do i=1,nx
+    do j=1,ny
+      do k=1,nz
+
+!!$direction represents the cylindrical radius here.
+        if (CCTK_EQUALS(shocktube_type,"xshock")) then
+           direction = sqrt((x(i,j,k)-shock_xpos)**2+&
+                (y(i,j,k)-shock_ypos)**2)
+        else if (CCTK_EQUALS(shocktube_type,"yshock")) then
+           direction = sqrt((y(i,j,k)-shock_ypos)**2+&
+                (z(i,j,k)-shock_zpos)**2)
+        else if (CCTK_EQUALS(shocktube_type,"zshock")) then
+           direction = sqrt((x(i,j,k)-shock_xpos)**2+&
+                (z(i,j,k)-shock_zpos)**2)
+        end if
+
+        tmp = cyl_fr(direction,r_inner,r_outer,rhol,rhor)
+        tmp2 = cyl_fr(direction,r_inner,r_outer,pressl,pressr)/( (gam - 1.d0) * tmp ) 
+        rho(i,j,k) = tmp
+        eps(i,j,k) = tmp2
+
+        velx(i,j,k) = 0.d0
+        vely(i,j,k) = 0.d0
+        velz(i,j,k) = 0.d0
+        Bvecx(i,j,k)=bvcxl
+        Bvecy(i,j,k)=bvcyl
+        Bvecz(i,j,k)=bvczl
+
+        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_CylindricalExplosionM
+
+
+
+function cyl_fr(R,r_inner,r_outer,min_f,max_f)
+
+  implicit none
+
+  CCTK_REAL :: cyl_fr
+  CCTK_REAL :: R,r_inner,r_outer,min_f,max_f
+
+  if(R .gt. r_outer) then 
+     cyl_fr = min_f
+  else if( R .lt. r_inner) then 
+     cyl_fr = max_f
+  else 
+     cyl_fr = exp( ( log(max_f)*(r_outer - R) + log(min_f)*(R - r_inner) ) / (r_outer - r_inner) ) 
+  end if
+
+  return 
+
+end function cyl_fr
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