merge branch hot_and_MHD_temp_dev into branch at revision r185

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

1 files changed, 137 insertions, 80 deletions

diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index e13c165..4e2631e 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -41,26 +41,49 @@ subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
 
   frac = CCTK_DELTA_TIME
 
-  !$OMP PARALLEL DO PRIVATE(j,i)
-  do k = 1, cctk_lsh(3)
-    do j = 1, cctk_lsh(2)
-      do i = 1, cctk_lsh(1)
-        if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. &
-             (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
-                                          type_bits, atmosphere)) .or. &
-             (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. &
-             (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
-          densrhs(i,j,k) = 0.0d0
-          srhs(i,j,k,:)   = 0.0d0
-          taurhs(i,j,k)  = 0.0d0
-          atmosphere_mask(i,j,k) = 1
-
-          SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
-        end if
-      end do
-    end do
-  end do
-  !$OMP END PARALLEL DO
+  if(evolve_temper.ne.1) then
+     !$OMP PARALLEL DO PRIVATE(j,i)
+     do k = 1, cctk_lsh(3)
+        do j = 1, cctk_lsh(2)
+           do i = 1, cctk_lsh(1)
+              if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. &
+                   (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
+                 type_bits, atmosphere)) .or. &
+                      (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. &
+                      (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
+                 densrhs(i,j,k) = 0.0d0
+                 srhs(i,j,k,:)   = 0.0d0
+                 taurhs(i,j,k)  = 0.0d0
+                 atmosphere_mask(i,j,k) = 1
+                 
+                 SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
+              end if
+           end do
+        end do
+     end do
+     !$OMP END PARALLEL DO
+  else
+     ! allow negative total energy
+     !$OMP PARALLEL DO PRIVATE(j,i)
+     do k = 1, cctk_lsh(3)
+        do j = 1, cctk_lsh(2)
+           do i = 1, cctk_lsh(1)
+              if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. &
+                   (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
+                 type_bits, atmosphere)) .or. &
+                      (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
+                 densrhs(i,j,k) = 0.0d0
+                 srhs(i,j,k,:)   = 0.0d0
+                 taurhs(i,j,k)  = 0.0d0
+                 atmosphere_mask(i,j,k) = 1
+                 
+                 SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
+              end if
+           end do
+        end do
+     end do
+     !$OMP END PARALLEL DO
+  endif
 
 end subroutine GRHydroUpdateAtmosphereMask
 
@@ -100,6 +123,10 @@ subroutine GRHydro_SetupMask(CCTK_ARGUMENTS)
      cons_tracerrhs = 0.0d0
   endif
 
+  if (evolve_y_e .ne. 0) then
+     y_e_con_rhs = 0.0d0
+  endif
+
   atmosphere_mask = 0
   
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
@@ -192,7 +219,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
   CCTK_REAL :: det, psi4pt
 
   CCTK_INT :: type_bits, atmosphere, not_atmosphere
-  
+
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
                               "in_atmosphere")
@@ -236,6 +263,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
     end do
   end do
 
+
 !!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
 
 end subroutine GRHydro_AtmosphereReset
@@ -254,14 +282,6 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
   CCTK_INT :: type_bits, atmosphere, not_atmosphere
   CCTK_INT :: eos_handle
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -272,7 +292,6 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   eos_handle = GRHydro_polytrope_handle
   
@@ -293,26 +312,40 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
         
         if (rho(i,j,k) .le. rho_min) then
           rho(i,j,k) = rho_min
+          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))
+
           velx(i,j,k) = 0.0d0
           vely(i,j,k) = 0.0d0
           velz(i,j,k) = 0.0d0
-#if USE_EOS_OMNI
-          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
-          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                 rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
-#else
-          press(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps(i,j,k))
-          eps(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, rho_min, press(i,j,k))
-#endif
-          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))
+
+          if(evolve_temper.ne.0) then
+             ! set the temperature to be relatively low
+             temperature(i,j,k) = 0.1d0
+             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))
+          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)
+             call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                  rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+             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))
+          endif
         end if
         if (timelevels .gt. 1) then
           if (rho_p(i,j,k) .le. rho_min) then
@@ -320,23 +353,38 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
             velx_p(i,j,k) = 0.0d0
             vely_p(i,j,k) = 0.0d0
             velz_p(i,j,k) = 0.0d0
-#if USE_EOS_OMNI
-          call EOS_Omni_press(eos_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)
-          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                 rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr)
-#else
-            press_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p(i,j,k))
-            eps_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, rho_min, press_p(i,j,k))
-#endif
+
             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))
+
+            if(evolve_temper.ne.0) then
+             ! set the temperature to be relatively low
+             temperature_p(i,j,k) = 0.1d0
+             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),&
+                  press_p(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_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),&
+                  temperature_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)
+               call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                    rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr)
+               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
+
           endif
         end if
         if (timelevels .gt. 2) then
@@ -345,23 +393,35 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
             velx_p_p(i,j,k) = 0.0d0
             vely_p_p(i,j,k) = 0.0d0
             velz_p_p(i,j,k) = 0.0d0
-#if USE_EOS_OMNI
-          call EOS_Omni_press(eos_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)
-          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                 rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr)
-#else
-            press_p_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p_p(i,j,k))
-            eps_p_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, rho_min, press_p_p(i,j,k))
-#endif
             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))
+
+            if(evolve_temper.ne.0) then
+             ! set the temperature to be relatively low
+             temperature_p_p(i,j,k) = 0.1d0
+             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),&
+                  press_p_p(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_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),&
+                  temperature_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)
+               call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                    rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr)
+               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
         endif
 
@@ -369,8 +429,5 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
     end do
   end do
 
-  write(*,*) "     GRHydro_InitialAtmosphereReset"
-!!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
-
 end subroutine GRHydro_InitialAtmosphereReset