[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 249)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Sat Jun 11 16:51:15 CDT 2011
User: cott
Date: 2011/06/11 04:51 PM
Modified:
/trunk/src/
GRHydro_Con2Prim.F90, GRHydro_Prim2Con.F90, GRHydro_UpdateMask.F90
Log:
* improve atmosphere handling so that it is now possible to use
hot EOS with an atmosphere
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +26 -11
===================================================================
--- trunk/src/GRHydro_Con2Prim.F90 2011-05-23 22:31:18 UTC (rev 248)
+++ trunk/src/GRHydro_Con2Prim.F90 2011-06-11 21:51:15 UTC (rev 249)
@@ -88,7 +88,7 @@
GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
!$OMP PARALLEL DO PRIVATE(i,j,itracer,&
- !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative)
+ !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp)
do k = 1, nz
do j = 1, ny
do i = 1, nx
@@ -128,7 +128,6 @@
Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
endif
-
if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
@@ -137,15 +136,26 @@
vel(i,j,k,:) = 0.d0
w_lorentz(i,j,k) = 1.d0
- 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)
+ if(evolve_temper.ne.0) then
+ ! set hot atmosphere values
+ temperature(i,j,k) = grhydro_hot_atmo_temp
+ y_e(i,j,k) = grhydro_hot_atmo_Y_e
+ y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k)
+ 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)
+ else
+ ! use polytropic EOS
+ 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),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+ call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+ endif
- ! w_lorentz=1, so the expression for tau reduces to:
tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
-
+
cycle
end if
@@ -196,6 +206,7 @@
endif
endif
+
if (epsnegative) then
#if 0
@@ -848,10 +859,14 @@
rho = GRHydro_rho_min
udens = rho
dens = sqrt(det) * rho
-! epsilon = (sqrt( (utau + pnew + udens)**2) - pnew - udens)/udens
- epsilon = epsminl
+ temp = GRHydro_hot_atmo_temp
+ ye = GRHydro_hot_atmo_Y_e
+ keytemp=1
+ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
+ rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+ keytemp=0
! w_lorentz=1, so the expression for utau reduces to:
- utau = rho + rho*epsminl - udens
+ utau = rho + rho*epsilon - udens
sx = 0.d0
sy = 0.d0
sz = 0.d0
File [modified]: GRHydro_Prim2Con.F90
Delta lines: +0 -1
===================================================================
--- trunk/src/GRHydro_Prim2Con.F90 2011-05-23 22:31:18 UTC (rev 248)
+++ trunk/src/GRHydro_Prim2Con.F90 2011-06-11 21:51:15 UTC (rev 249)
@@ -303,7 +303,6 @@
endif
endif
-
vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
File [modified]: GRHydro_UpdateMask.F90
Delta lines: +50 -23
===================================================================
--- trunk/src/GRHydro_UpdateMask.F90 2011-05-23 22:31:18 UTC (rev 248)
+++ trunk/src/GRHydro_UpdateMask.F90 2011-06-11 21:51:15 UTC (rev 249)
@@ -219,13 +219,19 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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)
More information about the Commits
mailing list