[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 332)

[roland.haas at physics.gatech.edu]

User: rhaas
Date: 2012/05/14 10:37 AM

Modified:
 /trunk/src/
  GRHydro_Reconstruct.F90

Log:
 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.

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_Reconstruct.F90
Delta lines: +33 -8
===================================================================
--- trunk/src/GRHydro_Reconstruct.F90	2012-05-14 15:31:15 UTC (rev 331)
+++ trunk/src/GRHydro_Reconstruct.F90	2012-05-14 15:37:35 UTC (rev 332)
@@ -44,6 +44,7 @@
 
   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 @@
     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,11 +126,16 @@
                    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
      if(evolve_mhd.ne.0) then
@@ -133,3 +157,4 @@
 
 end subroutine Reconstruction
 
+