diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/GRHydro_ShockTube.F90 | 153 |
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 |