GRHydro: Fixes to atmosphere mask and the way we sync.

 These are the changes as discussed with Roland. In
 particular, we sync the atmo mask in an extra call
 before all other syncs and then do atmo reset.
 Furthermore, primitives must always be synced in last
 post step in order to get initial guesses in the
 buffer zones via prolongation. Atmo reset must be
 done before prolongation to ensure we don't overwrite
 "good" data from the coarse grid!

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@377 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

5 files changed, 124 insertions, 6 deletions

diff --git a/interface.ccl b/interface.ccl
index 25f3d60..79ce24c 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -424,6 +424,14 @@ int GRHydro_atmosphere_mask type = GF Timelevels = 1 tags='Prolongation="None"'
   atmosphere_mask
 } "Flags to say whether a point needs to be reset to the atmosphere"
 
+# This real mask is set during UpdateAtmosphereMask and sync'ed afterwards (including possible interpatch interpolation)
+# After syn'ing and before any con2prim fun, we set the integer mask above based on the real-valued mask.
+# This ensures that any routine using the int mask is still correctly working.
+real GRHydro_atmosphere_mask_real type = GF Timelevels = 1 tags='Prolongation="sync" checkpoint="no"'
+{
+  atmosphere_mask_real
+} "Flags to say whether a point needs to be reset to the atmosphere. This is sync'ed (and possibly interpolated)!"
+
 int GRHydro_atmosphere_descriptors type=SCALAR
 {
   atmosphere_field_descriptor
diff --git a/schedule.ccl b/schedule.ccl
index 1bd3721..713a094 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -477,6 +477,7 @@ schedule GRHydro_Scalar_Setup IN CCTK_POSTINITIAL BEFORE HydroBase_RHS
 if (CCTK_EQUALS(evolution_method, "GRHydro")) {
 
 STORAGE:GRHydro_atmosphere_mask
+STORAGE:GRHydro_atmosphere_mask_real
 STORAGE:GRHydro_atmosphere_descriptors
 
 schedule GRHydro_SetupMask AT CCTK_Initial BEFORE HydroBase_Initial
@@ -999,6 +1000,25 @@ schedule group Do_GRHydro_Boundaries IN HydroBase_Boundaries
 {
 } "GRHydro Boundary conditions group"
 
+
+# Synchronize the atmosphere mask before doing other boundaries!
+schedule group GRHydro_AtmosphereMaskBoundaries IN HydroBase_PostStep BEFORE (HydroBase_Boundaries,GRHydro_PrimitiveBoundaries)
+{
+} "Apply boundary conditions to primitives"
+
+schedule GRHydro_SelectAtmosphereMaskBoundaries IN GRHydro_AtmosphereMaskBoundaries
+{
+  LANG: Fortran
+  OPTIONS: LEVEL
+  SYNC: GRHydro_atmosphere_mask_real
+} "Select primitive variables for boudary conditions"
+
+schedule group ApplyBCs AS GRHydro_ApplyAtmosphereMaskBCs in GRHydro_AtmosphereMaskBoundaries AFTER GRHydro_SelectAtmosphereMaskBoundaries
+{
+} "Apply boundary conditions to real-valued atmosphere mask"
+
+
+# Now synchronize other evolution variables
 if (!sync_conserved_only)
 {
    schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
@@ -1050,7 +1070,9 @@ else
     {
     } "Apply boundary conditions to primitives"
     
-    # This is necessary to provide initial guesses for Con2Prim after atmo reset (which happens only at last MoL post-step)!
+    # Primitive sync is always necessary to provide initial guesses for Con2Prim in the buffer zones via prolongation.
+    # Buffer zones are not guaranteed to have the correct values! Prolongation happens only in the last MoL poststep
+    # so we only need to sync in the last MoL post step.
     schedule group GRHydro_PrimitiveBoundaries IN HydroBase_PostStep BEFORE HydroBase_Boundaries IF GRHydro::InLastMoLPostStep
     {
     } "Apply boundary conditions to primitives"
@@ -1160,14 +1182,21 @@ schedule GRHydro_ClearLastMoLPostStep IN CCTK_WRAGH
 ### This is executed in the last MoL_PosStep                 ###
 ################################################################
 
+# First (right after sync) we need to convert real (and synchronized mask) to an integer mask!
+schedule GRHydroPostSyncAtmosphereMask IN HydroBase_PostStep AFTER GRHydro_AtmosphereMaskBoundaries
+{
+  LANG: Fortran
+} "Set integer atmosphere mask from synchronized real atmosphere mask"
+
+# Afterwards, reset mask!
 if (CCTK_Equals(Bvec_evolution_method,"GRHydro"))
 {
-  schedule GRHydro_AtmosphereResetM IN MoL_PostStep AFTER ADMBase_SetADMVars AFTER GRHydro_SetLastMoLPostStep BEFORE HydroBase_PostStep IF GRHydro::InLastMoLPostStep
+  schedule GRHydro_AtmosphereResetM IN HydroBase_PostStep AFTER GRHydroPostSyncAtmosphereMask BEFORE (HydroBase_Boundaries,GRHydro_PrimitiveBoundaries) IF GRHydro::InLastMoLPostStep
   {
     LANG: Fortran
   } "Reset the atmosphere - MHD version"
 } else {
-  schedule GRHydro_AtmosphereReset IN MoL_PostStep AFTER ADMBase_SetADMVars BEFORE HydroBase_PostStep IF GRHydro::InLastMoLPostStep
+  schedule GRHydro_AtmosphereReset IN HydroBase_PostStep AFTER GRHydroPostSyncAtmosphereMask BEFORE (HydroBase_Boundaries,GRHydro_PrimitiveBoundaries) IF GRHydro::InLastMoLPostStep
   {
     LANG: Fortran
   } "Reset the atmosphere"
diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90
index b3d153d..61f6664 100644
--- a/src/GRHydro_Boundaries.F90
+++ b/src/GRHydro_Boundaries.F90
@@ -99,6 +99,8 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
     !endif
   endif
   
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::atmosphere_mask_real")
+  
   sym(1) = -1
   sym(2) = 1
   sym(3) = 1
@@ -490,3 +492,49 @@ subroutine GRHydro_SelectPrimitiveBoundaries(CCTK_ARGUMENTS)
 
 end subroutine GRHydro_SelectPrimitiveBoundaries
 
+
+
+
+
+
+subroutine GRHydro_SelectAtmosphereMaskBoundaries(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer, dimension(3) :: sw
+  integer :: ierr = 0
+  integer :: i,j,k
+
+  CCTK_INT, parameter :: faces=CCTK_ALL_FACES
+  CCTK_INT, parameter :: ione=1
+  
+  sw = GRHydro_stencil
+
+
+!!$Flat boundaries if required  
+
+  if (verbose.eq.1) call CCTK_INFO("Selecting atmosphere mask BC")
+
+
+  if (CCTK_EQUALS(bound,"flat")) then
+      ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::GRHydro_atmosphere_mask_real", "Flat")
+  endif
+
+  if (CCTK_EQUALS(bound,"none")) then      
+      ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::GRHydro_atmosphere_mask_real", "None")
+  end if
+
+  if (CCTK_EQUALS(bound,"scalar")) then
+    call CCTK_WARN(0, "Until somebody uses this I see no reason to support it")
+  end if
+
+  if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition")
+
+end subroutine GRHydro_SelectAtmosphereMaskBoundaries
+
diff --git a/src/GRHydro_Eigenproblem.F90 b/src/GRHydro_Eigenproblem.F90
index 03576dd..fd9decc 100644
--- a/src/GRHydro_Eigenproblem.F90
+++ b/src/GRHydro_Eigenproblem.F90
@@ -88,7 +88,7 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps, &
   if(cs2.lt.0.0d0) then
      !$OMP CRITICAL
      if (abs(cs2) .gt. 1.0d-4) then
-        write(warnline,'(a50,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho
+        write(warnline,'(a50,6g16.7)') 'rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho
         call CCTK_WARN(1,warnline)
         call CCTK_WARN(1,"cs2 < 0! Check speed of sound calculation!")
         cs2 = 0.0d0
diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index 1808df3..19cacf3 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -48,7 +48,9 @@ subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
                  densrhs(i,j,k) = 0.0d0
                  srhs(i,j,k,:)   = 0.0d0
                  taurhs(i,j,k)  = 0.0d0
-                 atmosphere_mask(i,j,k) = 1
+                 ! Set real-valued mask! This will be sync'ed and right after syncing translated to
+                 ! our standard integer based mask (so that atmosphere_mask is still valid!).
+                 atmosphere_mask_real(i,j,k) = 1
               end if
            end do
         end do
@@ -67,7 +69,9 @@ subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
                  densrhs(i,j,k) = 0.0d0
                  srhs(i,j,k,:)   = 0.0d0
                  taurhs(i,j,k)  = 0.0d0
-                 atmosphere_mask(i,j,k) = 1
+                 ! Set real-valued mask! This will be sync'ed and right after syncing translated to
+                 ! our standard integer based mask (so that atmosphere_mask is still valid!).
+                 atmosphere_mask_real(i,j,k) = 1
               end if
            end do
         end do
@@ -77,6 +81,33 @@ subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
 
 end subroutine GRHydroUpdateAtmosphereMask
 
+
+subroutine GRHydroPostSyncAtmosphereMask(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: i,j,k
+
+!! This sets the integer atmo mask based on the real-valued (and sync'ed) atmo mask
+
+   !$OMP PARALLEL DO PRIVATE(k,j,i)
+   do k = 1, cctk_lsh(3)
+      do j = 1, cctk_lsh(2)
+	 do i = 1, cctk_lsh(1)
+	    if ( abs(atmosphere_mask_real(i,j,k)) .gt. 0.5) then
+	       atmosphere_mask(i,j,k) = 1
+	    end if
+	 end do
+      end do
+   end do
+   !$OMP END PARALLEL DO
+
+end subroutine GRHydroPostSyncAtmosphereMask
+
+
  /*@@
    @routine    GRHydro_SetupMask
    @date       Thu Jun 20 13:27:28 2002
@@ -116,6 +147,7 @@ subroutine GRHydro_SetupMask(CCTK_ARGUMENTS)
   endif
 
   atmosphere_mask = 0
+  atmosphere_mask_real = 0
   
   call CCTK_INFO("Setting up the atmosphere mask: all points are not_atmosphere")
   
@@ -263,6 +295,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
                   velz(i,j,k), 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
+                atmosphere_mask_real(i, j, k) = 0
              end if
           endif