add excision

however: be carefull with ENO and excision for the time being


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

8 files changed, 334 insertions, 285 deletions

diff --git a/src/Whisky_Con2Prim.F90 b/src/Whisky_Con2Prim.F90
index 5902c77..19d9a5c 100644
--- a/src/Whisky_Con2Prim.F90
+++ b/src/Whisky_Con2Prim.F90
@@ -146,14 +146,13 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   character(len=100) warnline
   
   CCTK_INT :: type_bits, atmosphere
-  CCTK_INT :: type2_bits, excised
+  CCTK_INT :: type2_bits
 
   CCTK_REAL :: local_min_tracer
 
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-    type2_bits = -1
-    excised = -1
+  type2_bits = -1
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -175,8 +174,9 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
     do j = 1, ny 
       do i = 1, nx
 
-        !do not compute if in atmosphere
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+        !do not compute if in atmosphere or in excised region
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+            (hydro_excision_mask(i,j,k) .ne. 0)) cycle
          
         epsnegative = .false.
 
@@ -481,6 +481,14 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
 
     f = pnew - EOS_Pressure(handle, rho, epsilon)
 
+    !check whether rho is NaN
+    if (rho .ne. rho) then      
+      !$OMP CRITICAL
+      write(warnline,'(a20,g20.7)') 'rho is NaN: ',rho
+      call CCTK_WARN(Whisky_NaN_verbose,warnline)
+      !$OMP END CRITICAL
+      goto 51
+    end if
     
   enddo
   
@@ -594,7 +602,7 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
   character(len=100) warnline
  
   CCTK_INT :: type_bits, atmosphere
-  CCTK_INT :: type2_bits, excised
+  CCTK_INT :: type2_bits
 
   CCTK_REAL :: local_min_tracer
   
@@ -608,8 +616,7 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
                                                                               
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-    type2_bits = -1
-    excised = -1
+  type2_bits = -1
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -625,8 +632,9 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
     do j = whisky_stencil, ny - whisky_stencil + 1
       do i = whisky_stencil, nx - whisky_stencil + 1
 
-        !do not compute if in atmosphere
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+        !do not compute if in atmosphere or in an excised region
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+            (hydro_excision_mask(i,j,k) .ne. 0)) cycle
          
         gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
         gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
@@ -809,7 +817,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt
 
   CCTK_INT :: type_bits, atmosphere
-  CCTK_INT :: type2_bits, excised
+  CCTK_INT :: type2_bits
 
   CCTK_REAL :: local_min_tracer
 !  character(len=400) :: warnline
@@ -817,8 +825,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
                                                                                
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-    type2_bits = -1
-    excised = -1
+  type2_bits = -1
  
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -839,9 +846,9 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
     do j = 1, ny
       do i = 1, nx
 
-        !do not compute if in atmosphere
-        !or in an excised region
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+        !do not compute if in atmosphere or in an excised region
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+            (hydro_excision_mask(i,j,k) .ne. 0)) cycle
 
         if (conformal_state .eq. 0) then
           call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
@@ -929,7 +936,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 #undef _EOS_BASE_INC_
 #endif
 #include "EOS_Base.inc"
-  
 
 !!$  Undensitize the variables 
 
@@ -954,7 +960,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     !$OMP CRITICAL
     write(warnline,'(a70,8g15.6)') 'NaN produced in sqrt(): (dens, det, s2, udens, enthalpy, x, y, z)', dens, det, s2, udens, enthalpy, x, y, z
     call CCTK_WARN(Whisky_NaN_verbose, warnline)
-    call CCTK_WARN(0, warnline)
     !$OMP END CRITICAL
   endif
 !!$ END: Check for NaN value (1st check)
@@ -1018,6 +1023,12 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     rhoold = rhonew
     rhonew = rhoold - f/df
 
+    if (rhonew .lt. 0.d0) then
+      whisky_C2P_failed = 1
+      !$OMP CRITICAL
+      call CCTK_WARN(Whisky_NaN_verbose, 'rhonew went below 0')
+      !$OMP END CRITICAL
+    end if
     
 !!$    Recalculate primitive variables and function
        
@@ -1145,12 +1156,11 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS)
        gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l
 
   CCTK_INT :: type_bits, atmosphere
-  CCTK_INT :: type2_bits, excised
+  CCTK_INT :: type2_bits
                                                                                
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-    type2_bits = -1
-    excised = -1
+  type2_bits = -1
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -1160,8 +1170,9 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS)
     do j = whisky_stencil, ny - whisky_stencil + 1
       do i = whisky_stencil, nx - whisky_stencil + 1
         
-        !do not compute if in atmosphere
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+        !do not compute if in atmosphere or in an excised region
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+            (hydro_excision_mask(i,j,k) .ne. 0)) cycle
         
         gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
         gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
@@ -1447,7 +1458,7 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
   logical, dimension(:,:,:), allocatable :: point_loop_condition
 
   CCTK_INT :: type_bits, atmosphere
-  CCTK_INT :: type2_bits, excised
+  CCTK_INT :: type2_bits
 
   CCTK_REAL :: local_min_tracer
 
@@ -1463,7 +1474,6 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
   type2_bits = -1
-  excised = -1
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -1482,11 +1492,6 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
 
 !!$  Set up rho and epsilon
 
-!!!  write(6,"(A15,i3,1P10E18.9)") "Xatmo 34 4 4: ", atmosphere_mask(34,4,4), rho(34,4,4), &
-!!!       whisky_rho_min, &
-!!!       dens(34,4,4), sqrt(Whisky_Det(34,4,4))*whisky_rho_min*(1.0d0+whisky_atmo_tolerance)
-!!!  call flush(6)
-
   do k = 1, nz
     do j = 1, ny
       do i = 1, nx
@@ -1518,11 +1523,7 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
 
            SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
            atmosphere_mask(i,j,k) = 1
-!           call Whisky_DumpPointInAtmos(rho(i,j,k), &
-!                vel(i,j,k,1), vel(i,j,k,2), vel(i,j,k,3), & 
-!                eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & 
-!                dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), tau(i,j,k), & 
-!                Whisky_det(i,j,k), whisky_rho_min)
+
         end if
         
         udens = dens(i,j,k) / sqrt(Whisky_det(i,j,k))
@@ -1652,13 +1653,13 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
       do j = 1, ny 
         do i = 1, nx
 
-          !do not compute if in atmosphere
-          !or in an excised region
-
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) &
+          !do not compute if in atmosphere or in an excised region
+          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+             (hydro_excision_mask(i,j,k) .ne. 0)) &
                then
              press_new(i,j,k) = press(i,j,k)
              press_old(i,j,k) = press(i,j,k)
+             cycle
           endif
 
 
@@ -1756,8 +1757,9 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
       do j = 1, ny 
         do i = 1, nx
 
-          !do not compute if in atmosphere
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+          !do not compute if in atmosphere or in an excised region
+          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+              (hydro_excision_mask(i,j,k) .ne. 0)) cycle
  
           f = press_new(i,j,k) - press(i,j,k)
 
@@ -1809,8 +1811,9 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
     do j = 1, ny 
       do i = 1, nx
         
-        !do not compute if in atmosphere
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+        !do not compute if in atmosphere or in an excised region
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+            (hydro_excision_mask(i,j,k) .ne. 0)) cycle
 
 !!$         Note that in the following we enforce eps > 0
 !!$         This makes the later check redundant
@@ -1903,11 +1906,6 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
 
   deallocate(f)
 
-!!!  write(6,"(A15,i3,1P10E18.9)") "Yatmo 34 4 4: ", atmosphere_mask(34,4,4), rho(34,4,4), &
-!!!       whisky_rho_min, &
-!!!       dens(34,4,4), sqrt(Whisky_Det(34,4,4))*whisky_rho_min*(1.0d0+whisky_atmo_tolerance)
-!!!  call flush(6)
-  
 end subroutine Conservative2PrimitiveGeneral
 
 
@@ -1949,14 +1947,13 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS)
   
   CCTK_INT :: type_bits
   CCTK_INT :: atmosphere
-  CCTK_INT :: type2_bits, excised
+  CCTK_INT :: type2_bits
 
   CCTK_REAL :: local_min_tracer
                                                                                
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-    type2_bits = -1
-    excised = -1
+  type2_bits = -1
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -2044,8 +2041,9 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS)
       do j = 1, ny 
         do i = 1, nx
 
-          !do not compute if in atmosphere
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+          !do not compute if in atmosphere or in an excised region
+          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+              (hydro_excision_mask(i,j,k) .ne. 0)) cycle
 
           enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
 
@@ -2096,8 +2094,9 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS)
     do j = 1, ny 
       do i = 1, nx
 
-        !do not compute if in atmosphere
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle
+        !do not compute if in atmosphere or in an excised region
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
+            (hydro_excision_mask(i,j,k) .ne. 0)) cycle
  
         udens = dens(i,j,k) / sqrt(Whisky_det(i,j,k))
 
diff --git a/src/Whisky_ENOReconstruct.F90 b/src/Whisky_ENOReconstruct.F90
index 093e524..42f0661 100644
--- a/src/Whisky_ENOReconstruct.F90
+++ b/src/Whisky_ENOReconstruct.F90
@@ -153,7 +153,7 @@ end subroutine Whisky_ENOShutdown
   (iand(mask((i)),(type_bits)).eq.(state_bits))
 
 subroutine Whisky_ENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
-     space_mask, excision_descriptors)
+     hydro_excision_mask)
   USE Whisky_ENOScalars
 
   implicit none
@@ -164,17 +164,13 @@ subroutine Whisky_ENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
   CCTK_REAL, dimension(nx) :: v, vplus, vminus
   CCTK_REAL, dimension(order, 1-order:nx+order) :: vdiff
 
-  CCTK_INT, dimension(nx) :: space_mask
+  CCTK_INT, dimension(nx) :: hydro_excision_mask
   logical, dimension(nx) :: trivial_rp
   logical, dimension(nx) :: excise
-  CCTK_INT :: excision_bits, excision_mask
-  CCTK_INT, dimension(3) :: excision_descriptors
+  logical :: normal_eno
 
   CCTK_REAL :: large = 1.d10
 
-  excision_bits=excision_descriptors(1)
-  excision_mask=excision_descriptors(2)
-
   vminus = 0.d0
   vplus = 0.d0
   vdiff = 0.d0
@@ -183,39 +179,81 @@ subroutine Whisky_ENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
   excise = .false.
   trivial_rp = .false.
 
-    
+!!$    Initialize excision
   do i = 1, nx
+    if (hydro_excision_mask(i) .ne. 0) then
+      vdiff(1, i) = large * (-1.d0*(1.d0+mod(i,10)))**i
+      trivial_rp(i) = .true.
+      excise(i) = .true.
+      if (i > 1) then
+        trivial_rp(i-1) = .true.
+      end if
+    end if
+  end do
 
+  do i = 1, nx
+!!$    Handle excision
+    normal_eno = .true.
+    if (i < nx) then
+     if (excise(i+1)) then
+      vminus(i) = v(i)
+      vplus(i) = v(i)
+      normal_eno = .false.
+     end if
+    end if
+    if (i > 1) then
+     if (excise(i-1)) then
+      vminus(i) = v(i)
+      vplus(i) = v(i)
+      normal_eno = .false.
+     end if
+    end if
+    if (normal_eno) then
 
 !!$    Calculate the undivided differences
-
       do k = 2, order
         do j = max(1, i - order), min(nx, i + order)
           vdiff(k, j) = vdiff(k - 1, j + 1) - vdiff(k - 1, j)
         end do
       end do
-      
+
 !!$    Ensure the stencil stays within the grid
-      
       vdiff(:, 1 - order : 0) = 1.d10
       vdiff(:, nx + 1 : nx + order) = 1.d10
-      
+
 !!$    Find the stencil
-      
       r = 0
-      
       do j = 2, order
         if ( abs(vdiff(j, i-r-1)) < abs(vdiff(j, i-r)) ) r = r + 1
       end do
       
 !!$    Calculate the reconstruction
-      
       do j = 0, order - 1
         vminus(i) = vminus(i) + eno_coeffs(r-1, j) * vdiff(1, i-r+j)
         vplus(i) = vplus(i) + eno_coeffs(r, j) * vdiff(1, i-r+j)
       end do
+    end if
+  end do
 
-    
+  do i = 1, nx
+    if (excise(i)) then
+      if (i > 1) then
+        if (.not. excise(i-1)) then
+          vminus(i) = vplus(i-1)
+        end if
+      end if
+      vplus(i) = vminus(i)
+    end if
+  end do
+  do i = nx, 1, -1
+    if (excise(i)) then
+      if (i < nx) then
+        if (.not. excise(i+1)) then
+          vplus(i) = vminus(i+1)
+        end if
+      end if
+      vminus(i) = vplus(i)
+    end if
   end do
 
 end subroutine Whisky_ENOReconstruct1d
diff --git a/src/Whisky_PPM.F90 b/src/Whisky_PPM.F90
index e249831..e56b910 100644
--- a/src/Whisky_PPM.F90
+++ b/src/Whisky_PPM.F90
@@ -53,7 +53,7 @@ end subroutine PPM_TVD
 subroutine SimplePPM_1d(handle,poly,&
      nx,dx,rho,velx,vely,velz,eps,press,rhominus,&
      velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,&
-     velzplus,epsplus,trivial_rp,space_mask, excision_descriptors,&
+     velzplus,epsplus,trivial_rp, hydro_excision_mask,&
      gxx, gxy, gxz, gyy, gyz, gzz, psi4, beta, alp, w_lorentz, &
      dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw)
 
@@ -89,9 +89,7 @@ subroutine SimplePPM_1d(handle,poly,&
 
   logical, dimension(nx) :: trivial_rp
 
-  CCTK_INT, dimension(nx) :: space_mask
-  CCTK_INT :: excision_bits, excision_mask
-  CCTK_INT, dimension(3) :: excision_descriptors
+  CCTK_INT, dimension(nx) :: hydro_excision_mask
 
   CCTK_REAL, dimension(nx) :: gxx, gxy, gxz, gyy, gyz, gzz, &
                               psi4, beta, alp, w_lorentz
@@ -432,7 +430,99 @@ do i = whisky_stencil, nx - whisky_stencil + 1
   end if
 
   !!$ excision
-  return
+ if (hydro_excision) then
+  do i = 1, nx
+    if (hydro_excision_mask(i) .ne. 0) then
+      if (i .gt. 1) then
+        trivial_rp(i-1)=.true.
+      end if
+      trivial_rp(i)=.true.
+    else
+      !!$ Do not optimize cond away by combining the 'if's. Fortran does not
+      !!$  have to follow the order of sub-expressions given here and might
+      !!$  access outside the array range
+      cond = i .gt. 1
+      if (cond) cond = hydro_excision_mask(i-1) .ne. 0
+      if (cond) then
+        rhominus(i)=rho(i)
+        rhoplus(i)=rho(i)
+        velxminus(i)=velx(i)
+        velxplus(i)=velx(i)
+        velyminus(i)=vely(i)
+        velyplus(i)=vely(i)
+        velzminus(i)=velz(i)
+        velzplus(i)=velz(i)
+        rhominus(i-1)=rho(i)
+        rhoplus(i-1)=rho(i)
+        velxminus(i-1)=velx(i)
+        velxplus(i-1)=velx(i)
+        velyminus(i-1)=vely(i)
+        velyplus(i-1)=vely(i)
+        velzminus(i-1)=velz(i)
+        velzplus(i-1)=velz(i)
+        if (poly .eq. 0) then
+          epsminus(i)=eps(i)
+          epsplus(i)=eps(i)
+          epsminus(i-1)=eps(i)
+          epsplus(i-1)=eps(i)
+        end if
+      else
+        cond = (i.gt.2) .and. (i.lt.nx)
+        if (cond) cond = (ppm_mppm .eq. 0) .and.&
+             (hydro_excision_mask(i-2) .ne. 0)
+        if (cond) then
+          call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
+          call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i))
+          call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i))
+          call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i))
+          if (poly .eq. 0) then
+            call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
+          end if
+        end if
+      end if
+      cond = i.lt.nx
+      if (cond) cond = hydro_excision_mask(i+1) .ne. 0;
+      if (cond) then
+        rhominus(i)=rho(i)
+        rhoplus(i)=rho(i)
+        velxminus(i)=velx(i)
+        velxplus(i)=velx(i)
+        velyminus(i)=vely(i)
+        velyplus(i)=vely(i)
+        velzminus(i)=velz(i)
+        velzplus(i)=velz(i)
+        rhominus(i+1)=rho(i)
+        rhoplus(i+1)=rho(i)
+        velxminus(i+1)=velx(i)
+        velxplus(i+1)=velx(i)
+        velyminus(i+1)=vely(i)
+        velyplus(i+1)=vely(i)
+        velzminus(i+1)=velz(i)
+        velzplus(i+1)=velz(i)
+        if (poly .eq. 0) then
+          epsminus(i)=eps(i)
+          epsplus(i)=eps(i)
+          epsminus(i+1)=eps(i)
+          epsplus(i+1)=eps(i)
+        endif
+      else
+        cond = (i.lt.nx-1) .and. (i.gt.1)
+        if (cond) cond = (ppm_mppm .eq. 0) .and.&
+             (hydro_excision_mask(i+2) .ne. 0)
+        if (cond) then
+          call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
+          call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i))
+          call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i))
+          call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i))
+          if (poly .eq. 0) then
+            call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
+          end if
+        end if
+      end if
+    end if
+  end do
+ end if
+return
 
 end subroutine SimplePPM_1d
 
diff --git a/src/Whisky_Prim2Con.F90 b/src/Whisky_Prim2Con.F90
index 01f3d9f..2c823f5 100644
--- a/src/Whisky_Prim2Con.F90
+++ b/src/Whisky_Prim2Con.F90
@@ -329,10 +329,13 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
   
   implicit none
   
+  DECLARE_CCTK_PARAMETERS
+
   CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
   CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
        deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet   
   CCTK_INT :: handle
+  character(len=256) NaN_WarnLine
 
 #ifdef _EOS_BASE_INC_
 #undef _EOS_BASE_INC_
@@ -348,6 +351,10 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
     ! and hard to trace wrong values below. There is no good value to
     ! choose in this case, but something small is probably the best of
     ! all bad choices.
+    !$OMP CRITICAL
+    write(NaN_WarnLine,'(a100,1g15.6)') 'Infinite Lorentz factor reset. rho: ', drho
+    call CCTK_WARN(Whisky_NaN_verbose, NaN_WarnLine)
+    !$OMP END CRITICAL
     w = 1.d-20
   else
     w = 1.d0 / sqrt(1.d0 - w_tmp)
diff --git a/src/Whisky_Reconstruct.F90 b/src/Whisky_Reconstruct.F90
index 89bedb7..7b2365d 100644
--- a/src/Whisky_Reconstruct.F90
+++ b/src/Whisky_Reconstruct.F90
@@ -55,7 +55,6 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
        &type_bitsy, trivialy, not_trivialy, &
        &type_bitsz, trivialz, not_trivialz
 
-  CCTK_INT, dimension(3) :: excision_descriptors
   CCTK_REAL, dimension(:,:,:),allocatable :: &
        &psi4, lbetax, lbetay, lbetaz
 !!$  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: &
@@ -78,10 +77,6 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
     call CCTK_WARN(0, "Allocation problems with lbeta")
   end if
   
-    excision_descriptors(1)=-1
-    excision_descriptors(2)=-1
-    excision_descriptors(3)=-1
-
   call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX")
   call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", &
        &"trivial")
@@ -147,42 +142,32 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
            call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
                 tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
                 tracerminus(:,:,:,itracer), &
-                trivial_rp, space_mask, excision_descriptors)
+                trivial_rp, hydro_excision_mask)
         enddo
      end if
 
     if (CCTK_EQUALS(recon_vars,"primitive")) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           rho, rhoplus, rhominus, trivial_rp, space_mask, &
-           excision_descriptors)
+           rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           vel(:,:,:,1), velxplus, velxminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           vel(:,:,:,2), velyplus, velyminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           vel(:,:,:,3), velzplus, velzminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           eps, epsplus, epsminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           eps, epsplus, epsminus, trivial_rp, hydro_excision_mask)
     else if (CCTK_EQUALS(recon_vars,"conservative")) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           dens, densplus, densminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           dens, densplus, densminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           scon(:,:,:,1), sxplus, sxminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           scon(:,:,:,2), syplus, syminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           scon(:,:,:,3), szplus, szminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           tau, tauplus, tauminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           tau, tauplus, tauminus, trivial_rp, hydro_excision_mask)
     else
       call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
     end if
@@ -224,8 +209,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
                press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),&
                velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),&
                velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),&
-               trivial_rp(:,j,k),space_mask(:,j,k), &
-               excision_descriptors,&
+               trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
                gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
                gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),&
                w_lorentz(:,j,k), &
@@ -264,8 +248,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
                press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),&
                velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),&
                velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),&
-               trivial_rp(j,:,k),space_mask(j,:,k), &
-               excision_descriptors,&
+               trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
                gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
                gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),&
                w_lorentz(j,:,k), &
@@ -304,8 +287,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
                press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),&
                velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),&
                velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),&
-               trivial_rp(j,k,:),space_mask(j,k,:), &
-               excision_descriptors,&
+               trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
                gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
                gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),&
                w_lorentz(j,k,:), &
@@ -346,7 +328,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
            call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
                 tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
                 tracerminus(:,:,:,itracer), trivial_rp, &
-                space_mask, excision_descriptors)
+                hydro_excision_mask)
         enddo
      end if
 
@@ -357,45 +339,35 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
           if (CCTK_EQUALS(recon_vars,"primitive")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  dens(:,j,k),densminus(:,j,k),densplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  sy(:,j,k),syminus(:,j,k),syplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  sz(:,j,k),szminus(:,j,k),szplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
@@ -418,45 +390,35 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
           if (CCTK_EQUALS(recon_vars,"primitive")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  dens(j,:,k),densminus(j,:,k),densplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  sy(j,:,k),syminus(j,:,k),syplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  sz(j,:,k),szminus(j,:,k),szplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
@@ -479,45 +441,35 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
           if (CCTK_EQUALS(recon_vars,"primitive")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  dens(j,k,:),densminus(j,k,:),densplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  sy(j,k,:),syminus(j,k,:),syplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  sz(j,k,:),szminus(j,k,:),szplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
diff --git a/src/Whisky_ReconstructPoly.F90 b/src/Whisky_ReconstructPoly.F90
index dc065ce..c7a5dfa 100644
--- a/src/Whisky_ReconstructPoly.F90
+++ b/src/Whisky_ReconstructPoly.F90
@@ -55,7 +55,6 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
        &type_bitsy, trivialy, not_trivialy, &
        &type_bitsz, trivialz, not_trivialz
 
-  CCTK_INT, dimension(3) :: excision_descriptors
   CCTK_REAL, dimension(:,:,:),allocatable :: &
        &psi4, lbetax, lbetay, lbetaz
 !!$  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: &
@@ -78,10 +77,6 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
     call CCTK_WARN(0, "Allocation problems with lbeta")
   end if
   
-    excision_descriptors(1)=-1
-    excision_descriptors(2)=-1
-    excision_descriptors(3)=-1
-
   call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX")
   call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", &
        &"trivial")
@@ -147,38 +142,30 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
            call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
                 tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
                 tracerminus(:,:,:,itracer), &
-                trivial_rp, space_mask, excision_descriptors)
+                trivial_rp, hydro_excision_mask)
         enddo
      end if
      
 
     if (CCTK_EQUALS(recon_vars,"primitive")) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           rho, rhoplus, rhominus, trivial_rp, space_mask, &
-           excision_descriptors)
+           rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           vel(:,:,:,1), velxplus, velxminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           vel(:,:,:,2), velyplus, velyminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           vel(:,:,:,3), velzplus, velzminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
 
     else if (CCTK_EQUALS(recon_vars,"conservative")) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           dens, densplus, densminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           dens, densplus, densminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           scon(:,:,:,1), sxplus, sxminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           scon(:,:,:,2), syplus, syminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           scon(:,:,:,3), szplus, szminus, trivial_rp, space_mask, &
-           excision_descriptors)
+           scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask)
 
     else
       call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
@@ -221,8 +208,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
                press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),&
                velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),&
                velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),&
-               trivial_rp(:,j,k),space_mask(:,j,k), &
-               excision_descriptors,&
+               trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
                gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
                gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),&
                w_lorentz(:,j,k), &
@@ -261,8 +247,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
                press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),&
                velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),&
                velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),&
-               trivial_rp(j,:,k),space_mask(j,:,k), &
-               excision_descriptors,&
+               trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
                gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
                gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),&
                w_lorentz(j,:,k), &
@@ -301,8 +286,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
                press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),&
                velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),&
                velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),&
-               trivial_rp(j,k,:),space_mask(j,k,:), &
-               excision_descriptors,&
+               trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
                gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
                gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),&
                w_lorentz(j,k,:), &
@@ -343,7 +327,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
            call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
                 tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
                 tracerminus(:,:,:,itracer), trivial_rp, &
-                space_mask, excision_descriptors)
+                hydro_excision_mask)
         enddo
      end if
 
@@ -354,37 +338,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
           if (CCTK_EQUALS(recon_vars,"primitive")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  dens(:,j,k),densminus(:,j,k),densplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  sy(:,j,k),syminus(:,j,k),syplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  sz(:,j,k),szminus(:,j,k),szplus(:,j,k),&
-                 trivial_rp(:,j,k), space_mask(:,j,k), &
-                 excision_descriptors)
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
@@ -407,37 +383,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
           if (CCTK_EQUALS(recon_vars,"primitive")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  dens(j,:,k),densminus(j,:,k),densplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  sy(j,:,k),syminus(j,:,k),syplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  sz(j,:,k),szminus(j,:,k),szplus(j,:,k),&
-                 trivial_rp(j,:,k), space_mask(j,:,k), &
-                 excision_descriptors)
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
@@ -460,37 +428,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
           if (CCTK_EQUALS(recon_vars,"primitive")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  dens(j,k,:),densminus(j,k,:),densplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  sy(j,k,:),syminus(j,k,:),syplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
             call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  sz(j,k,:),szminus(j,k,:),szplus(j,k,:),&
-                 trivial_rp(j,k,:), space_mask(j,k,:), &
-                 excision_descriptors)
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
diff --git a/src/Whisky_TVDReconstruct.F90 b/src/Whisky_TVDReconstruct.F90
index d6f4b95..9788a05 100644
--- a/src/Whisky_TVDReconstruct.F90
+++ b/src/Whisky_TVDReconstruct.F90
@@ -14,8 +14,6 @@
 
 #include "SpaceMask.h"
 
-/* #define WHISKY_UNSTABLE_TVD */
-
  /*@@
    @routine    tvdreconstruct
    @date       Sat Jan 26 02:12:12 2002
@@ -32,7 +30,7 @@
 @@*/
 
 subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-     orig, bextp, bextm, trivial_rp, space_mask, excision_descriptors)
+     orig, bextp, bextm, trivial_rp, hydro_excision_mask)
   
   USE Whisky_Scalars
   
@@ -45,9 +43,7 @@ subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
   CCTK_REAL, dimension(nx, ny, nz) :: orig, bextp, bextm
   CCTK_REAL :: dupw, dloc, delta, ratio, hdelta
   logical, dimension(nx,ny,nz) :: trivial_rp
-  CCTK_INT, dimension(nx,ny,nz) :: space_mask
-  CCTK_INT :: excision_bits, excision_mask
-  CCTK_INT, dimension(3) :: excision_descriptors
+  CCTK_INT, dimension(nx,ny,nz) :: hydro_excision_mask
 
 
   bextp = 0.d0
@@ -60,41 +56,45 @@ subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
     do k = whisky_stencil, nz-whisky_stencil+1
       do j = whisky_stencil, ny-whisky_stencil+1
         do i = whisky_stencil, nx-whisky_stencil+1
-          dupw = orig(i, j, k) - orig(i-xoffset, j-yoffset, k-zoffset)
-          dloc = orig(i+xoffset, j+yoffset, k+zoffset) - orig(i, j, k)
-!!$        For minmod set everything here.
-!!$        Otherwise call the limiter function.
-        
-          if (MINMOD) then
-!!$            if (dupw*dloc < 0.d0) then
-!!$              delta=0.d0
-!!$            else if (abs(dupw)<abs(dloc)) then
-!!$              delta=dupw
-!!$            else
-!!$              delta=dloc
-!!$            end if
-
-            delta = minmod_func(dupw,dloc)
-
-!!$        This is an alternative equivalent implementation 
-!!$          of vanLeeer MC slopelimiter
-          else if (MC2) then
-            if (dupw*dloc < 0.d0) then
-              delta=0.d0
-            else 
-              delta=sign(min(2.d0*abs(dupw),2.d0*abs(dloc),&
-                   0.5d0*(abs(dupw)+abs(dloc))),dupw+dloc)
+          if (hydro_excision_mask(i,j,k) .ne. 0) then
+            trivial_rp(i-xoffset, j-yoffset, k-zoffset) = .true.
+            trivial_rp(i, j, k) = .true.
+            bextm(i, j, k) = orig(i, j, k)
+            bextp(i, j, k) = orig(i, j, k)
+            if (hydro_excision_mask(i+xoffset,j+yoffset,k+zoffset) .eq. 0) then
+              bextm(i, j, k) = orig(i+xoffset, j+yoffset, k+zoffset)
+              bextp(i, j, k) = orig(i+xoffset, j+yoffset, k+zoffset)
             end if
+          else if ((hydro_excision_mask(i-xoffset,j-yoffset,k-zoffset) .ne. 0) .or. &
+                   (hydro_excision_mask(i+xoffset,j+yoffset,k+zoffset) .ne. 0)) then
+            bextm(i, j, k) = orig(i, j, k)
+            bextp(i, j, k) = orig(i, j, k)
           else
-            delta = 0.5d0*(dupw + dloc)
-            if (abs(dupw) < myfloor ) dupw = myfloor*sign(1.d0, dupw) 
-            if (abs(dloc) < myfloor ) dloc = myfloor*sign(1.d0, dloc) 
-            ratio = dupw / dloc
-            call slopelimiter(ratio, delta)
+            dupw = orig(i, j, k) - orig(i-xoffset, j-yoffset, k-zoffset)
+            dloc = orig(i+xoffset, j+yoffset, k+zoffset) - orig(i, j, k)
+
+            if (MINMOD) then
+              delta = minmod_func(dupw,dloc)
+            else if (MC2) then
+!!$            This is an alternative equivalent implementation 
+!!$              of vanLeeer MC slopelimiter
+              if (dupw*dloc < 0.d0) then
+                delta=0.d0
+              else 
+                delta=sign(min(2.d0*abs(dupw),2.d0*abs(dloc),&
+                     0.5d0*(abs(dupw)+abs(dloc))),dupw+dloc)
+              end if
+            else
+              delta = 0.5d0*(dupw + dloc)
+              if (abs(dupw) < myfloor ) dupw = myfloor*sign(1.d0, dupw) 
+              if (abs(dloc) < myfloor ) dloc = myfloor*sign(1.d0, dloc) 
+              ratio = dupw / dloc
+              call slopelimiter(ratio, delta)
+            end if
+            hdelta = 0.5d0 * delta 
+            bextm(i, j, k) = orig(i, j, k) - hdelta
+            bextp(i, j, k) = orig(i, j, k) + hdelta
           end if
-          hdelta = 0.5d0 * delta 
-          bextm(i, j, k) = orig(i, j, k) - hdelta
-          bextp(i, j, k) = orig(i, j, k) + hdelta
         end do
       end do
     end do
diff --git a/src/Whisky_UpdateMask.F90 b/src/Whisky_UpdateMask.F90
index 49e28c3..0b590ff 100644
--- a/src/Whisky_UpdateMask.F90
+++ b/src/Whisky_UpdateMask.F90
@@ -34,10 +34,12 @@ subroutine WhiskyUpdateAtmosphereMask(CCTK_ARGUMENTS)
 
   frac = CCTK_DELTA_TIME
 
+  !$OMP PARALLEL DO PRIVATE(j,i)
   do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
       do i = 1, cctk_lsh(1)
-        if ( (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
+        if ( (hydro_excision_mask(i,j,k) .ne. 0) .or. &
+             (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
                                           type_bits, atmosphere)) .or. &
              (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. &
              (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
@@ -51,6 +53,7 @@ subroutine WhiskyUpdateAtmosphereMask(CCTK_ARGUMENTS)
       end do
     end do
   end do
+  !$OMP END PARALLEL DO
 
 end subroutine WhiskyUpdateAtmosphereMask