diff options
Diffstat (limited to 'src/GRHydro_UpdateMaskM.F90')
| -rw-r--r-- | src/GRHydro_UpdateMaskM.F90 | 205 |
1 files changed, 155 insertions, 50 deletions
diff --git a/src/GRHydro_UpdateMaskM.F90 b/src/GRHydro_UpdateMaskM.F90 index d6a6bb8..8c3feef 100644 --- a/src/GRHydro_UpdateMaskM.F90 +++ b/src/GRHydro_UpdateMaskM.F90 @@ -44,11 +44,18 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt +! 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 + do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - - if ( atmosphere_mask(i, j, k) .eq. 1) then + + if (atmosphere_mask(i, j, k) .ne. 0) then rho(i,j,k) = GRHydro_rho_min vel(i,j,k,1) = 0.0d0 @@ -56,20 +63,43 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) vel(i,j,k,3) = 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 prim2conpolytypeM(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), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& - rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), & - vel(i,j,k,3), eps(i,j,k), press(i,j,k), & - Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k)) - if (wk_atmosphere .eq. 0) then - atmosphere_mask(i, j, k) = 0 - atmosphere_mask_real(i, j, k) = 0 - 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 prim2conM_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),Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),& + eps(i,j,k),press(i,j,k),Bvec(i,j,k,1), & + Bvec(i,j,k,2), Bvec(i,j,k,3), 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 prim2conpolytypeM(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), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), & + vel(i,j,k,3), eps(i,j,k), press(i,j,k), & + Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k)) + if (wk_atmosphere .eq. 0) then + atmosphere_mask(i, j, k) = 0 + end if + + end if + endif end do end do @@ -88,6 +118,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt + CCTK_REAL :: rho_min CCTK_INT :: eos_handle @@ -100,23 +131,50 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) eos_handle = GRHydro_polytrope_handle + rho_min = GRHydro_rho_min + if (initial_atmosphere_factor .gt. 0) then + rho_min = rho_min * initial_atmosphere_factor + endif + 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 + if (rho(i,j,k) .le. rho_min) then + rho(i,j,k) = rho_min vel(i,j,k,1) = 0.0d0 vel(i,j,k,2) = 0.0d0 vel(i,j,k,3) = 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)) + + 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 prim2conM_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),Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),& + eps(i,j,k),press(i,j,k),Bvec(i,j,k,1), & + Bvec(i,j,k,2), Bvec(i,j,k,3), 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(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) - 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 prim2conpolytypeM(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, & @@ -125,6 +183,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), & vel(i,j,k,3), eps(i,j,k), press(i,j,k), & Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k)) + end if end if if (timelevels .gt. 1) then if (rho_p(i,j,k) .le. GRHydro_rho_min) then @@ -132,47 +191,93 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) vel_p(i,j,k,1) = 0.0d0 vel_p(i,j,k,2) = 0.0d0 vel_p(i,j,k,3) = 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)) + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + 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),& + press_p(i,j,k),keyerr,anyerr) + call prim2conM_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),Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& + rho_p(i,j,k),vel_p(i,j,k,1),vel_p(i,j,k,2),vel_p(i,j,k,3),& + eps_p(i,j,k),press_p(i,j,k),Bvec_p(i,j,k,1), & + Bvec_p(i,j,k,2), Bvec_p(i,j,k,3), 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) - 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 - 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 prim2conpolytypeM(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), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& - rho_p(i,j,k), vel_p(i,j,k,1), vel_p(i,j,k,2), & - vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), & - Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),w_lorentz_p(i,j,k)) - endif + 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) + call prim2conpolytypeM(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), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& + rho_p(i,j,k), vel_p(i,j,k,1), vel_p(i,j,k,2), & + vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), & + Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),w_lorentz_p(i,j,k)) + + end if + end if 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 vel_p_p(i,j,k,1) = 0.0d0 vel_p_p(i,j,k,2) = 0.0d0 vel_p_p(i,j,k,3) = 0.0d0 - 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) - 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 prim2conpolytypeM(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), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& - rho_p_p(i,j,k), vel_p_p(i,j,k,1), vel_p_p(i,j,k,2), & - vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), & - Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),w_lorentz_p_p(i,j,k)) - endif - endif + gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k)) + + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + 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),& + press_p_p(i,j,k),keyerr,anyerr) + call prim2conM_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),Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& + rho_p_p(i,j,k),vel_p_p(i,j,k,1),vel_p_p(i,j,k,2),vel_p_p(i,j,k,3),& + eps_p_p(i,j,k),press_p_p(i,j,k),Bvec_p_p(i,j,k,1), & + Bvec_p_p(i,j,k,2), Bvec_p_p(i,j,k,3), 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(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) + call prim2conpolytypeM(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), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& + rho_p_p(i,j,k), vel_p_p(i,j,k,1), vel_p_p(i,j,k,2), & + vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), & + Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),w_lorentz_p_p(i,j,k)) + + end if + end if + end if end do end do |