Update grid.cc

grid.cc

@@ -2533,12 +2533,18 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
            (grid::get_cellcoordmin(cellindex, d) - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) -
           tstart;
       d_coordminboundary[d] = CLIGHT_PROP * t_coordminboundary;
+      if (last_cross == posdirections[d]) {
+        d_coordminboundary[d] = -1;
+      }
 
       const double t_coordmaxboundary = ((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) /
                                          (cellcoordmax[d] - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) -
                                         tstart;
 
       d_coordmaxboundary[d] = CLIGHT_PROP * t_coordmaxboundary;
+      if (last_cross == negdirections[d]) {
+        d_coordmaxboundary[d] = -1;
+      }
     }
   } else {
     assert_always(false);
@@ -2550,7 +2556,7 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
   int snext = 0;
   for (int d = 0; d < get_ndim(GRID_TYPE); d++) {
     // upper d coordinate of the current cell
-    if ((d_coordmaxboundary[d] > 0) && (d_coordmaxboundary[d] < distance) && (last_cross != negdirections[d])) {
+    if ((d_coordmaxboundary[d] > 0) && (d_coordmaxboundary[d] < distance)) {
       choice = posdirections[d];
       distance = d_coordmaxboundary[d];
       if (grid::get_cellcoordpointnum(cellindex, d) == (grid::ncoordgrid[d] - 1)) {
@@ -2562,7 +2568,7 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
     }
 
     // lower d coordinate of the current cell
-    if ((d_coordminboundary[d] > 0) && (d_coordminboundary[d] < distance) && (last_cross != posdirections[d])) {
+    if ((d_coordminboundary[d] > 0) && (d_coordminboundary[d] < distance)) {
       choice = negdirections[d];
       distance = d_coordminboundary[d];
       if (grid::get_cellcoordpointnum(cellindex, d) == 0) {