Snapshot of ET GRHydro_InitData rev. 132 with a poloidal magnetic field routine

added.

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

1 files changed, 166 insertions, 0 deletions

diff --git a/src/GRHydro_RotorM.F90 b/src/GRHydro_RotorM.F90
new file mode 100644
index 0000000..6191b8a
--- /dev/null
+++ b/src/GRHydro_RotorM.F90
@@ -0,0 +1,166 @@
+ /*@@
+   @file      GRHydro_RotorM.F90
+   @date      Aug 15, 2011
+   @author    Scott Noble, Joshua Faber, Bruno Mundim
+   @desc 
+   Cylindrical magnetized rotor test (see Etienne et al.).
+   @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)
+#define Bconsx(i,j,k) Bcons(i,j,k,1)
+#define Bconsy(i,j,k) Bcons(i,j,k,2)
+#define Bconsz(i,j,k) Bcons(i,j,k,3)
+
+
+ /*@@
+   @routine    GRHydro_Rotor
+   @date       Aug 11, 2011
+   @author     Scott Noble, Joshua Faber, Bruno Mundim
+   @desc 
+   Initial data for magnetic rotor - parameters from Etienne et al. 2010
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Using GRHydro_ShockTubeM.F90 as a template.
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_RotorM(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  CCTK_INT :: i, j, k, nx, ny, nz
+  CCTK_REAL :: radius, det, v_max
+  CCTK_REAL :: rhoin, rhoout, pressin, pressout
+  CCTK_REAL :: bvcxl,bvcyl,bvczl
+  CCTK_REAL :: tmp,tmp2,gam,r_rot
+
+  CCTK_INT :: use_smoothing
+  CCTK_REAL :: rsmooth_rel,rfact
+
+!!$Hardwired bfield and gamma values
+  bvcxl = 1.0
+  bvcyl = 0.0
+  bvczl = 0.0
+  
+!!$Adiabatic index for test
+  gam = (5.d0/3.d0)
+  
+!!$radius of rotor is 0.1
+  r_rot = 0.1
+  
+!!$Maximum velocity is 0.995
+  v_max = 0.995
+
+!!$Inner values
+  rhoin   = 10.d0
+  pressin = 1.d0
+  
+!!$Outer values
+  rhoout   = 1.d0
+  pressout = 1.d0
+
+
+!!$ Smooth the edge?  Define the radius in relative terms if so
+  use_smoothing = 0
+  rsmooth_rel = 0.05
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  if(v_max>1.0) then
+     write(*,*)"No superluminal speeds!!!!  v_max=",v_max," reset to 0.995"
+     v_max=0.995
+  endif
+
+  do i=1,nx
+     do j=1,ny
+        do k=1,nz
+           
+           radius = sqrt(x(i,j,k)**2+y(i,j,k)**2)
+           
+           if(radius.le.r_rot) then
+              rho(i,j,k) = rhoin
+              press(i,j,k) = pressin
+              velx(i,j,k) = v_max/r_rot*y(i,j,k)
+              vely(i,j,k) = -1.d0*v_max/r_rot*x(i,j,k)
+              velz(i,j,k) = 0.d0
+           else if((use_smoothing.eq.1) .and. &
+                ((radius.gt.r_rot) .and. &
+                (radius.le.((1.0+rsmooth_rel)*r_rot)))) then
+              rfact = (radius/r_rot - 1.d0) / rsmooth_rel
+              rho(i,j,k) = rfact*rhoout + (1.d0-rfact)*rhoin
+              press(i,j,k) = rfact*pressout + (1.d0-rfact)*pressin
+              velx(i,j,k) = (1.d0-rfact)*v_max * y(i,j,k) / radius
+              velx(i,j,k) = -1.d0*(1.d0-rfact)*v_max * x(i,j,k) / radius
+              velz(i,j,k) = 0.d0
+           else
+              rho(i,j,k) = rhoout
+              press(i,j,k) = pressout
+              velx(i,j,k) = 0.d0
+              vely(i,j,k) = 0.d0
+              velz(i,j,k) = 0.d0
+           endif
+           eps(i,j,k)=press(i,j,k)/(gam-1.d0)/rho(i,j,k)
+           
+           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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
+                   w_lorentz(i,j,k))
+           end if
+           
+        enddo
+     enddo
+  enddo
+  
+  densrhs = 0.d0
+  srhs = 0.d0
+  taurhs = 0.d0
+  Bconsrhs = 0.d0
+
+
+  return
+  
+end subroutine GRHydro_RotorM
+
+