* improve atmosphere handling so that it is now possible to use

  hot EOS with an atmosphere



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

1 files changed, 50 insertions, 23 deletions

diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index a097446..3f3be88 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -219,13 +219,19 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
   CCTK_REAL :: det, psi4pt
 
   CCTK_INT :: type_bits, atmosphere, not_atmosphere
+! begin EOS Omni vars
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress,xeps,xtemp,xye
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0
+! end EOS Omni vars                             
 
   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")
-
+!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr)
   do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
       do i = 1, cctk_lsh(1)
@@ -235,37 +241,52 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
                     atmosphere)) &
              &) then
 
-!!$          write(*,*) 'Resetting at ',i,j,k, atmosphere_mask(i, j, k), &
-!!$             & (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\
-!!$                    atmosphere)) 
-
           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(GRHydro_polytrope_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))
-          if (wk_atmosphere .eq. 0) then
-            atmosphere_mask(i, j, k) = 0
-            SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
-                                      not_atmosphere)
-          end if
 
-        end if
+          if(evolve_temper.ne.0) then
+!             ! set the temperature to be relatively low
+             temperature(i,j,k) = grhydro_hot_atmo_temp
+             y_e(i,j,k) = grhydro_hot_atmo_Y_e
+             keytemp = 1
+             call EOS_Omni_press(GRHydro_eos_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)
+
+             call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
+                  i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),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),&
+                  temperature(i,j,k),y_e(i,j,k))
+             y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k)
+          else
+             call prim2conpolytype(GRHydro_polytrope_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))
+             if (wk_atmosphere .eq. 0) then
+                atmosphere_mask(i, j, k) = 0
+                SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
+                not_atmosphere)
+             end if
+          endif
+
+       end if
 
       end do
     end do
   end do
+!$OMP END PARALLEL DO
 
 
-!!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
-
 end subroutine GRHydro_AtmosphereReset
 
 subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
@@ -316,8 +337,9 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
           velz(i,j,k) = 0.0d0
 
           if(evolve_temper.ne.0) then
-             ! set the temperature to be relatively low
-             temperature(i,j,k) = 0.1d0
+!             ! set the temperature to be relatively low
+             temperature(i,j,k) = grhydro_hot_atmo_temp
+             y_e(i,j,k) = grhydro_hot_atmo_Y_e
              keytemp = 1
              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
@@ -330,6 +352,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
                   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),&
                   temperature(i,j,k),y_e(i,j,k))
+             y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k)
           else
              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
@@ -355,7 +378,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
 
             if(evolve_temper.ne.0) then
              ! set the temperature to be relatively low
-             temperature_p(i,j,k) = 0.1d0
+             temperature_p(i,j,k) = grhydro_hot_atmo_temp
+             y_e_p(i,j,k) = grhydro_hot_atmo_Y_e
              keytemp = 1
              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho_p(i,j,k),eps_p(i,j,k),temperature_p(i,j,k),y_e_p(i,j,k),&
@@ -368,6 +392,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
                   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),&
                   temperature_p(i,j,k),y_e_p(i,j,k))
+             y_e_con_p(i,j,k) = dens_p(i,j,k) * y_e_p(i,j,k)
             else
                call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                     rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr)
@@ -394,7 +419,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
 
             if(evolve_temper.ne.0) then
              ! set the temperature to be relatively low
-             temperature_p_p(i,j,k) = 0.1d0
+             temperature_p_p(i,j,k) = grhydro_hot_atmo_temp
+             y_e_p_p(i,j,k) = grhydro_hot_atmo_Y_e
              keytemp = 1
              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),&
@@ -406,6 +432,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
                   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),&
                   temperature_p_p(i,j,k),y_e_p_p(i,j,k))
+             y_e_con_p_p(i,j,k) = dens_p_p(i,j,k) * y_e_p_p(i,j,k)
             else
                call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                     rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr)