* improve handling of cs2 < 0 in the "standard" (non-MHD, non-hot EOS) eigenvalue routine.

Patch by Christian Ott.


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

4 files changed, 61 insertions, 10 deletions

diff --git a/interface.ccl b/interface.ccl
index 25f3d60..689f155 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -424,6 +424,15 @@ 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="copy"'
+{
+  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..37cb55b 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
@@ -1024,6 +1025,7 @@ if (!sync_conserved_only)
       SYNC: Y_e_con
       SYNC: lvel
       SYNC: lBvec
+      SYNC: GRHydro_atmosphere_mask_real  # Must be synced otherwise we get inconsistent results when using multiple procs
    } "Select GRHydro boundary conditions"
 }
 else
@@ -1039,6 +1041,7 @@ else
       SYNC: psidc
       SYNC: GRHydro_cons_tracers
       SYNC: Y_e_con
+      SYNC: GRHydro_atmosphere_mask_real  # Must be synced otherwise we get inconsistent results when using multiple procs
     } "Select GRHydro boundary conditions"
 
     # after a regrid Cactus relies on the boundary thorns to fill in the outer
@@ -1050,11 +1053,6 @@ 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)!
-    schedule group GRHydro_PrimitiveBoundaries IN HydroBase_PostStep BEFORE HydroBase_Boundaries IF GRHydro::InLastMoLPostStep
-    {
-    } "Apply boundary conditions to primitives"
-    
 
     if(CCTK_IsImplementationActive("Coordinates")) {
     
@@ -1160,14 +1158,20 @@ schedule GRHydro_ClearLastMoLPostStep IN CCTK_WRAGH
 ### This is executed in the last MoL_PosStep                 ###
 ################################################################
 
+# First we need to convert real (and synchronized mask) to an integer mask!
+schedule GRHydroPostSyncAtmosphereMask IN HydroBase_PostStep AFTER HydroBase_Boundaries
+{
+  LANG: Fortran
+} "Set integer atmosphere mask from synchronized real atmosphere 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 BEFORE HydroBase_Con2Prim 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 BEFORE HydroBase_Con2Prim IF GRHydro::InLastMoLPostStep
   {
     LANG: Fortran
   } "Reset the atmosphere"
diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90
index b3d153d..227a00a 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
@@ -210,6 +212,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
          "GRHydro::tau", "Flat")
     ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
          "GRHydro::scon", "Flat")
+    ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::GRHydro_atmosphere_mask_real", "Flat")
     if (sync_conserved_only .eq. 0) then
       ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
 	    "HydroBase::w_lorentz", "Flat")
@@ -282,7 +286,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
          "GRHydro::tau", "None")
     ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
          "GRHydro::scon", "None")
-    
+    ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::GRHydro_atmosphere_mask_real", "None") 
     if (sync_conserved_only .eq. 0) then
       ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
 	    "HydroBase::w_lorentz", "None")
diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index 1808df3..c11d218 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 ( atmosphere_mask_real(i,j,k) .ne. 0) 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