* add a hot shocktube to run shocktube tests with the nuclear EOS

From: Christian Ott <cott@bethe.tapir.caltech.edu>

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

1 files changed, 153 insertions, 0 deletions

diff --git a/src/GRHydro_ShockTube.F90 b/src/GRHydro_ShockTube.F90
index 8b9105b..7ad5b54 100644
--- a/src/GRHydro_ShockTube.F90
+++ b/src/GRHydro_ShockTube.F90
@@ -164,3 +164,156 @@ subroutine GRHydro_shocktube(CCTK_ARGUMENTS)
   return
   
 end subroutine GRHydro_shocktube
+
+subroutine GRHydro_shocktube_hot(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, templ, tempr, yel, yer
+  CCTK_REAL :: vlowx, vlowy, vlowz, w, tenthalpy
+
+! begin EOS Omni vars                                                                                                  
+  integer :: n,keytemp,anyerr,keyerr(1)
+  integer :: handle
+  character(len=256) :: warnline
+  ! handle for nuclear EOS
+  handle=4
+  n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0
+! end EOS Omni vars      
+  
+  nx = cctk_ash(1)
+  ny = cctk_ash(2)
+  nz = cctk_ash(3)
+  
+  call CCTK_INFO("Setting up initial data for hot shocktube")
+
+  if(.not.CCTK_EQUALS(Y_e_evolution_method,"GRHydro").or.&
+       .not.CCTK_EQUALS(temperature_evolution_method,"GRHydro")) then
+     call CCTK_WARN(0,"Must have Y_e_evolution_method and temperature_evolution_method set to GRHydro")
+  endif
+
+  if (nuceos_read_table.eq.0) then
+     call CCTK_WARN(0,"You must read in a nuclear EOS table for initial data shocktube_hot to work!")
+  endif
+
+  if (.not.CCTK_EQUALS(GRHydro_eos_table,"nuc_eos").or..not.CCTK_EQUALS(GRHydro_eos_type,"General")) then
+     call CCTK_WARN(0,"You must set GRHydro::GRHydro_eos_table = nuc_eos and GRHydro::GRHydro_eos_type = General!")
+  endif
+
+  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 = 1.62e-8
+          rhor = 1.62e-9
+          velxl = 0.d0
+          velxr = 0.d0
+          velyl = 0.d0
+          velyr = 0.d0
+          velzl = 0.d0
+          velzr = 0.d0
+          templ = 8.0d0
+          tempr = 0.6d0
+          yel = 0.48d0
+          yer = 0.48d0
+        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
+          rho(i,j,k) = rhol
+          velx(i,j,k) = velxl
+          vely(i,j,k) = velyl
+          velz(i,j,k) = velzl
+          temperature(i,j,k) = templ
+          y_e(i,j,k) = yel
+        else
+          rho(i,j,k) = rhor
+          velx(i,j,k) = velxr
+          vely(i,j,k) = velyr
+          velz(i,j,k) = velzr
+          temperature(i,j,k) = tempr
+          y_e(i,j,k) = yer
+        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))
+
+        ! call EOS with
+        keytemp = 1
+        call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+             rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),press(i,j,k),keyerr,anyerr)
+
+        if(anyerr.ne.0) then
+           call CCTK_WARN(1,"Error in Initial Data EOS call!")
+           write(warnline,"(A10,i8)") "keyerr= ",keyerr
+           call CCTK_WARN(0,warnline)
+        endif
+
+        ! set up conserved variables
+        w = 1.0d0 / &
+             sqrt(1.0d0 - (gxx(i,j,k)*vel(i,j,k,1)*vel(i,j,k,1) &
+                + gyy(i,j,k)*vel(i,j,k,2)*vel(i,j,k,2) &
+                + gzz(i,j,k)*vel(i,j,k,3)*vel(i,j,k,3) ) )
+
+        vlowx = gxx(i,j,k)*vel(i,j,k,1) &
+             + gxy(i,j,k)*vel(i,j,k,2)  &
+             + gxz(i,j,k)*vel(i,j,k,3)
+        vlowy = gxy(i,j,k)*vel(i,j,k,1) &
+             + gyy(i,j,k)*vel(i,j,k,2)  &
+             + gyz(i,j,k)*vel(i,j,k,3)
+        vlowz = gxz(i,j,k)*vel(i,j,k,1) &
+             + gyz(i,j,k)*vel(i,j,k,2)  &
+             + gzz(i,j,k)*vel(i,j,k,3)
+
+
+        dens(i,j,k) = sqrt(det)*w*rho(i,j,k)
+
+        tenthalpy = 1.0d0 + eps(i,j,k) + press(i,j,k) / rho(i,j,k)
+
+        tau(i,j,k) = sqrt(det)*( (rho(i,j,k)*(1.0d0+eps(i,j,k))+press(i,j,k))*w*w - press(i,j,k)) &
+             - dens(i,j,k)
+
+        w_lorentz(i,j,k) = w
+
+        scon(i,j,k,1) = sqrt(det)*rho(i,j,k)*tenthalpy*(w**2) &
+             *vlowx
+        scon(i,j,k,2) = sqrt(det)*rho(i,j,k)*tenthalpy*(w**2) &
+             *vlowy
+        scon(i,j,k,3) = sqrt(det)*rho(i,j,k)*tenthalpy*(w**2) &
+             *vlowz
+        
+        Y_e_con(i,j,k) = dens(i,j,k)*Y_e(i,j,k)
+           
+     enddo
+    enddo
+  enddo
+
+  densrhs = 0.d0
+  srhs = 0.d0
+  taurhs = 0.d0
+
+  return
+  
+end subroutine GRHydro_shocktube_hot