From a8d0fcb877d375f4544460cec1bf8be6bccb30fa Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 14 May 2012 15:37:35 +0000 Subject: GRHydro: correct atmosphere handling after primitve reconstruction I have fixed a bug in GRHydro's atmosphere handling after reconstruction. After reconstruction, the reconstructed plus and minus face values are tested for whether they drop below atmosphere level and are then reset. In particular, this means that plus and minus face values for cell i don't generally coincide with the minus/plus face values of the corresponding cells i-1 and i+1. If the reconstruction previously found that the Riemann problem was trivial between, say rhominus(i+1) and rhoplus(i), when changing rhoplus(i), it is not anymore! However, the mask "trivial_rp" testing for a trivial Riemann problem is not changed! This leads to inconsistent behavior at the boundary of the atmosphere. The symptom was a drift in the center of mass of a static and perfectly symmetric TOV star. Patch by Christian Reisswig. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@332 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Reconstruct.F90 | 41 +++++++++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 8 deletions(-) diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index c35afec..f10bc71 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -44,6 +44,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) CCTK_INT :: i,j,k CCTK_REAL :: local_min_tracer + CCTK_INT :: type_bits, not_trivial ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -83,12 +84,30 @@ subroutine Reconstruction(CCTK_ARGUMENTS) end if end if + + if (flux_direction == 1) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(not_trivial, "Hydro_RiemannProblemX", & + &"not_trivial") + else if (flux_direction == 2) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(not_trivial, "Hydro_RiemannProblemY", & + &"not_trivial") + else if (flux_direction == 3) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(not_trivial, "Hydro_RiemannProblemZ", & + &"not_trivial") + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) - if(rhoplus(i,j,k).lt.GRHydro_rho_min .or. & - rhominus(i,j,k).lt.GRHydro_rho_min) then + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + if(rhoplus(i,j,k).lt.GRHydro_rho_min .or. & + rhominus(i,j,k).lt.GRHydro_rho_min) then + rhoplus(i,j,k) = rho(i,j,k) rhominus(i,j,k) = rho(i,j,k) epsplus(i,j,k) = eps(i,j,k) @@ -107,10 +126,15 @@ subroutine Reconstruction(CCTK_ARGUMENTS) tracerminus(i,j,k,:) = tracer(i,j,k,:) end where end if - enddo - enddo + ! Riemann problem might not be trivial anymore!!!! + SpaceMask_SetStateBitsF90(space_mask, i-xoffset, j-yoffset, k-zoffset, type_bits, not_trivial) + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_trivial) + enddo enddo - !$OMP END PARALLEL DO + enddo + !$OMP END PARALLEL DO + + if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then @@ -133,3 +157,4 @@ subroutine Reconstruction(CCTK_ARGUMENTS) end subroutine Reconstruction + -- cgit v1.2.3