diff options
Diffstat (limited to 'src/GRHydro_UpdateMaskM.F90')
| -rw-r--r-- | src/GRHydro_UpdateMaskM.F90 | 227 |
1 files changed, 227 insertions, 0 deletions
diff --git a/src/GRHydro_UpdateMaskM.F90 b/src/GRHydro_UpdateMaskM.F90 new file mode 100644 index 0000000..5e46380 --- /dev/null +++ b/src/GRHydro_UpdateMaskM.F90 @@ -0,0 +1,227 @@ + /*@@ + @file GRHydro_UpdateMaskM.F90 + @date Sep 2, 2010 + @author + @desc + Alter the update terms if inside the atmosphere or excision region + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +#include "GRHydro_Macros.h" +#include "SpaceMask.h" + +!!$ We don't need to adapt GRHydroUpdateAtmosphereMask, GRHydro_SetupMask, or +!!$ since we need to evolve Bvec in the atmosphere + +!!$ In GRHydro_AtmosphereResetM, we just need to switch the P2C calls to MHD + + /*@@ + @routine GRHydro_AtmosphereResetM + @date Sep 2, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + After MoL has evolved, if a point is supposed to be reset then do so. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i, j, k + 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") + call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& + "not_in_atmosphere") + + 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) & + &.or. (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\ + 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 + 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)) + 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), Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(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), 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 + + end do + end do + end do + +!!$ call GRHydro_BoundariesM(CCTK_PASS_FTOF) + +end subroutine GRHydro_AtmosphereResetM + +subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i, j, k + CCTK_REAL :: det, psi4pt + + 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 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 +! end EOS Omni vars +#endif + + eos_handle = GRHydro_polytrope_handle + + 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") + + 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 + vel(i,j,k,1) = 0.0d0 + vel(i,j,k,2) = 0.0d0 + vel(i,j,k,3) = 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, GRHydro_rho_min, eps(i,j,k)) + eps(i,j,k) = EOS_SpecificIntEnergy(eos_handle, GRHydro_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 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, & + dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + tau(i,j,k), Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(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), w_lorentz(i,j,k)) + end if + if (timelevels .gt. 1) then + if (rho_p(i,j,k) .le. GRHydro_rho_min) then + rho_p(i,j,k) = GRHydro_rho_min + vel_p(i,j,k,1) = 0.0d0 + vel_p(i,j,k,2) = 0.0d0 + vel_p(i,j,k,3) = 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, GRHydro_rho_min, eps_p(i,j,k)) + eps_p(i,j,k) = EOS_SpecificIntEnergy(eos_handle, GRHydro_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 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), Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_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), w_lorentz_p(i,j,k)) + endif + 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 +#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, GRHydro_rho_min, eps_p_p(i,j,k)) + eps_p_p(i,j,k) = EOS_SpecificIntEnergy(eos_handle, GRHydro_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 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), Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_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), w_lorentz_p_p(i,j,k)) + endif + endif + + end do + end do + end do + + write(*,*) " GRHydro_InitialAtmosphereReset" +!!$ call GRHydro_BoundariesM(CCTK_PASS_FTOF) + +end subroutine GRHydro_InitialAtmosphereResetM + |