let GRHydro set the atmosphere of initial data, so that ID thorns don't have to worry about it. This changes the data slitely for EOSG testsuites because this now uses EOSG calls instead of the otherwise fixed calculations of the TOVSolver, thus the change of two testsuites

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@156 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 91 insertions, 0 deletions

diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index d7cdff0..45dfce0 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -17,6 +17,12 @@
 #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 velx_p(i,j,k) vel_p(i,j,k,1)
+#define vely_p(i,j,k) vel_p(i,j,k,2)
+#define velz_p(i,j,k) vel_p(i,j,k,3)
+#define velx_p_p(i,j,k) vel_p_p(i,j,k,1)
+#define vely_p_p(i,j,k) vel_p_p(i,j,k,2)
+#define velz_p_p(i,j,k) vel_p_p(i,j,k,3)
 
 subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
 
@@ -233,3 +239,88 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
 !!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
 
 end subroutine GRHydro_AtmosphereReset
+
+subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: i, j, k
+  CCTK_REAL :: det, psi4pt
+
+  CCTK_INT :: type_bits, atmosphere, not_atmosphere
+  CCTK_INT :: eos_handle
+
+  eos_handle = GRHydro_polytrope_handle
+  if (use_eosgeneral == 0) then
+    eos_handle = -1
+  endif
+  
+  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
+  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
+                              "in_atmosphere")
+  call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
+                              "not_in_atmosphere")
+
+  do k = 1, cctk_lsh(3)
+    do j = 1, cctk_lsh(2)
+      do i = 1, cctk_lsh(1)
+        
+        if (rho(i,j,k) .le. GRHydro_rho_min) then
+          rho(i,j,k) = GRHydro_rho_min
+          velx(i,j,k) = 0.0d0
+          vely(i,j,k) = 0.0d0
+          velz(i,j,k) = 0.0d0
+          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 prim2conpolytype(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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
+               tau(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
+        if (timelevels .gt. 1) then
+          if (rho_p(i,j,k) .le. GRHydro_rho_min) then
+            rho_p(i,j,k) = GRHydro_rho_min
+            velx_p(i,j,k) = 0.0d0
+            vely_p(i,j,k) = 0.0d0
+            velz_p(i,j,k) = 0.0d0
+            det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \
+                  gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k))
+            call prim2conpolytype(eos_handle, &
+                 gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), &
+                 gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, &
+                 dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), &
+                 tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), &
+                 velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k))
+          endif
+        end if
+        if (timelevels .gt. 2) then
+          if (rho_p_p(i,j,k) .le. GRHydro_rho_min) then
+            rho_p_p(i,j,k) = GRHydro_rho_min
+            velx_p_p(i,j,k) = 0.0d0
+            vely_p_p(i,j,k) = 0.0d0
+            velz_p_p(i,j,k) = 0.0d0
+            det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \
+                  gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k))
+            call prim2conpolytype(eos_handle, &
+                 gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), &
+                 gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, &
+                 dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), &
+                 tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &
+                 velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k))
+          endif
+        endif
+
+      end do
+    end do
+  end do
+
+  write(*,*) "     GRHydro_InitialAtmosphereReset"
+!!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
+
+end subroutine GRHydro_InitialAtmosphereReset
+