From c141ed811d4a5b84ecc077999d1fec0a3c598b95 Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 14 Jan 2013 14:23:27 +0000 Subject: GRHydro: Introduce variable atmopshere level / tolerance as function of radius. Off by default. From: Christian Reisswig git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@453 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_UpdateMask.F90 | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) (limited to 'src/GRHydro_UpdateMask.F90') diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90 index 67a42cb..a93872b 100644 --- a/src/GRHydro_UpdateMask.F90 +++ b/src/GRHydro_UpdateMask.F90 @@ -198,14 +198,17 @@ subroutine GRHydro_InitAtmosMask(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS CCTK_INT :: i,j,k - - !$OMP PARALLEL DO PRIVATE(i,j,k) + CCTK_REAL :: dummy1, dummy2 + + !$OMP PARALLEL DO PRIVATE(i,j,k, dummy1,dummy2) 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 + !if (rho(i,j,k) .le. GRHydro_rho_min) then + IF_BELOW_ATMO(rho(i,j,k), GRHydro_rho_min, 0.0, r(i,j,k)) then atmosphere_mask(i,j,k) = 1 atmosphere_mask_real(i,j,k) = 1 end if @@ -247,7 +250,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: det, dummy1, dummy2 ! save memory when MP is not used @@ -283,14 +286,14 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) if (verbose.eq.1) call CCTK_INFO("Entering AtmosphereReset.") -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) if (atmosphere_mask(i, j, k) .ne. 0) then - rho(i,j,k) = GRHydro_rho_min + SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.0d0 @@ -359,7 +362,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det - CCTK_REAL :: rho_min + CCTK_REAL :: rho_min, dummy1, dummy2 CCTK_INT :: eos_handle @@ -447,15 +450,16 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) rho_min = rho_min * initial_atmosphere_factor endif -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if (rho(i,j,k) .le. rho_min .or. & + SET_ATMO_MIN(dummy2, rho_min, r(i,j,k)) + if (rho(i,j,k) .le. dummy2 .or. & GRHydro_enable_internal_excision /= 0 .and. & hydro_excision_mask(i,j,k) .ne. 0) then - rho(i,j,k) = rho_min + SET_ATMO_MIN(rho(i,j,k), dummy2, r(i,j,k)) det = SPATIAL_DETERMINANT(g11(i,j,k), g12(i,j,k), g13(i,j,k), \ g22(i,j,k), g23(i,j,k), g33(i,j,k)) @@ -496,8 +500,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) endif end if if (timelevels .gt. 1) then - if (rho_p(i,j,k) .le. rho_min) then - rho_p(i,j,k) = rho_min + if (rho_p(i,j,k) .le. dummy2) then + SET_ATMO_MIN(rho_p(i,j,k), dummy2, r(i,j,k)) velx_p(i,j,k) = 0.0d0 vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0 @@ -540,8 +544,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) endif end if if (timelevels .gt. 2) then - if (rho_p_p(i,j,k) .le. rho_min) then - rho_p_p(i,j,k) = rho_min + if (rho_p_p(i,j,k) .le. dummy2) then + SET_ATMO_MIN(rho_p_p(i,j,k), dummy2, r(i,j,k)) velx_p_p(i,j,k) = 0.0d0 vely_p_p(i,j,k) = 0.0d0 velz_p_p(i,j,k) = 0.0d0 -- cgit v1.2.3