diff --git a/grid.cc b/grid.cc
index 4ba333ce3..44ae9c452 100644
--- a/grid.cc
+++ b/grid.cc
@@ -2440,22 +2440,23 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
 
   const auto negdirections = std::array<enum cell_boundary, 3>{COORD0_MIN, COORD1_MIN, COORD2_MIN};
   const auto posdirections = std::array<enum cell_boundary, 3>{COORD0_MAX, COORD1_MAX, COORD2_MAX};
-  if constexpr (TESTMODE) {
+  if constexpr (TESTMODE || true) {
     for (int d = 0; d < get_ndim(GRID_TYPE); d++) {
       // flip is either zero or one to indicate +ve and -ve boundaries along the selected axis
       for (int flip = 0; flip < 2; flip++) {
         bool isoutside_thisside = false;
         double delta = 0.;
+        const bool pos_component_vel_relative_to_flow = (pktvelgridcoord[d] * tstart) > pktposgridcoord[d];
         if (flip != 0) {
           // packet pos below min
           const double boundaryposmin = grid::get_cellcoordmin(cellindex, d) / globals::tmin * tstart;
           delta = pktposgridcoord[d] - boundaryposmin;
-          isoutside_thisside = pktposgridcoord[d] < (boundaryposmin * 0.99999);  // accuracy tolerance
+          isoutside_thisside = pos_component_vel_relative_to_flow && (pktposgridcoord[d] < (boundaryposmin - 10.));
         } else {
           // packet pos above max
           const double boundaryposmax = cellcoordmax[d] / globals::tmin * tstart;
           delta = pktposgridcoord[d] - boundaryposmax;
-          isoutside_thisside = pktposgridcoord[d] > (boundaryposmax * 1.00001);
+          isoutside_thisside = !pos_component_vel_relative_to_flow && (pktposgridcoord[d] > (boundaryposmax + 10.));
         }
 
         if (isoutside_thisside) {