diff options
author | cott <cott@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-06-11 21:51:15 +0000 |
---|---|---|
committer | cott <cott@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-06-11 21:51:15 +0000 |
commit | 44a25b645da3bf416016be364f086d95e77c7ee8 (patch) | |
tree | c9621e08d7843dfb9a8c56fb7df60edf019d1f21 | |
parent | 81835fa64f41e7cd4e190158e591269c09fa829a (diff) |
* 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
-rw-r--r-- | src/GRHydro_Con2Prim.F90 | 39 | ||||
-rw-r--r-- | src/GRHydro_Prim2Con.F90 | 1 | ||||
-rw-r--r-- | src/GRHydro_UpdateMask.F90 | 73 |
3 files changed, 77 insertions, 36 deletions
diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 9a3d2bc..9d69ebd 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -88,7 +88,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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) - - 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) + 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) + 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 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) endif endif + if (epsnegative) then #if 0 @@ -848,10 +859,14 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve 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 diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 17548a0..0bfea69 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -303,7 +303,6 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & endif endif - vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz 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) |