Move checks into boundary_distance()

artis-mcrt · Mar 7, 2025 · 2a2296e · 2a2296e
1 parent 095cb37
commit 2a2296e
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 59 deletions.
diff --git a/gammapkt.cc b/gammapkt.cc
@@ -729,32 +729,7 @@ void transport_gamma(Packet &pkt, const double t2) {
   // boundaries. sdist is the boundary distance and snext is the
   // grid cell into which we pass.
 
-  auto [sdist, snext] = grid::boundary_distance(pkt.dir, pkt.pos, pkt.prop_time, pkt.where, &pkt.last_cross);
-
-  const double maxsdist = (GRID_TYPE == GridType::CARTESIAN3D)
-                              ? globals::rmax * pkt.prop_time / globals::tmin
-                              : 2 * globals::rmax * (pkt.prop_time + sdist / CLIGHT_PROP) / globals::tmin;
-  if (sdist > maxsdist) {
-    printout("Unreasonably large sdist (gamma). Abort. %g %g %g\n", globals::rmax, pkt.prop_time / globals::tmin,
-             sdist);
-    assert_always(false);
-  }
-
-  if (sdist < 0) {
-    printout("Negative distance (sdist). Abort?\n");
-    sdist = 0;
-  }
-
-  if (((snext < 0) && (snext != -99)) || (snext >= grid::ngrid)) {
-    printout("Heading for inappropriate grid cell. Abort.\n");
-    printout("Current cell %d, target cell %d.\n", pkt.where, snext);
-    assert_always(false);
-  }
-
-  if (sdist > globals::max_path_step) {
-    sdist = globals::max_path_step;
-    snext = pkt.where;
-  }
+  const auto [sdist, snext] = grid::boundary_distance(pkt.dir, pkt.pos, pkt.prop_time, pkt.where, &pkt.last_cross);
 
   // Now consider the scattering/destruction processes.
   // Compton scattering - need to determine the scattering co-efficient.

diff --git a/grid.cc b/grid.cc
@@ -4,6 +4,7 @@
 
 #include <algorithm>
 #include <array>
+#include <cassert>
 #include <cctype>
 #include <cmath>
 #include <cstddef>
@@ -2598,6 +2599,19 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
     }
   }
 
+  const double maxsdist = (GRID_TYPE == GridType::CARTESIAN3D)
+                              ? globals::rmax * tstart / globals::tmin
+                              : 2 * globals::rmax * (tstart + distance / CLIGHT_PROP) / globals::tmin;
+  assert_always(distance <= maxsdist);
+  assert_always(distance >= 0);
+
+  assert_always((snext == -99) || ((snext >= 0) && (snext < grid::ngrid)));
+
+  if (distance > globals::max_path_step) {
+    distance = globals::max_path_step;
+    snext = cellindex;
+  }
+
   return {distance, snext};
 }
 

diff --git a/rpkt.cc b/rpkt.cc
@@ -626,46 +626,14 @@ auto do_rpkt_step(Packet &pkt, const double t2) -> bool {
   // Start by finding the distance to the crossing of the grid cell
   // boundaries. sdist is the boundary distance and snext is the
   // grid cell into which we pass.
-  auto [sdist, snext] = grid::boundary_distance(pkt.dir, pkt.pos, pkt.prop_time, pkt.where, &pkt.last_cross);
+  const auto [sdist, snext] = grid::boundary_distance(pkt.dir, pkt.pos, pkt.prop_time, pkt.where, &pkt.last_cross);
 
   if (sdist == 0) {
     grid::change_cell(pkt, snext);
     const int new_nonemptymgi = grid::get_propcell_nonemptymgi(pkt.where);
 
     return (pkt.type == TYPE_RPKT && (new_nonemptymgi < 0 || new_nonemptymgi == nonemptymgi));
   }
-  const double maxsdist = (GRID_TYPE == GridType::CARTESIAN3D)
-                              ? globals::rmax * pkt.prop_time / globals::tmin
-                              : 2 * globals::rmax * (pkt.prop_time + sdist / CLIGHT_PROP) / globals::tmin;
-  if (sdist > maxsdist) {
-    printout("[fatal] do_rpkt: Unreasonably large sdist for packet %d. Rpkt. Abort. %g %g %g\n", pkt.number,
-             globals::rmax, pkt.prop_time / globals::tmin, sdist);
-    std::abort();
-  }
-
-  if (sdist < 0) {
-    const int cellindexnew = pkt.where;
-    printout("[warning] r_pkt: Negative distance (sdist = %g). Abort.\n", sdist);
-    printout("[warning] r_pkt: cell %d snext %d\n", cellindexnew, snext);
-    printout("[warning] r_pkt: pos %g %g %g\n", pkt.pos[0], pkt.pos[1], pkt.pos[2]);
-    printout("[warning] r_pkt: dir %g %g %g\n", pkt.dir[0], pkt.dir[1], pkt.dir[2]);
-    printout("[warning] r_pkt: cell corner %g %g %g\n",
-             grid::get_cellcoordmin(cellindexnew, 0) * pkt.prop_time / globals::tmin,
-             grid::get_cellcoordmin(cellindexnew, 1) * pkt.prop_time / globals::tmin,
-             grid::get_cellcoordmin(cellindexnew, 2) * pkt.prop_time / globals::tmin);
-    printout("[warning] r_pkt: cell width %g\n", grid::wid_init(cellindexnew, 0) * pkt.prop_time / globals::tmin);
-    assert_always(false);
-  }
-  if (((snext != -99) && (snext < 0)) || (snext >= grid::ngrid)) {
-    printout("[fatal] r_pkt: Heading for inappropriate grid cell. Abort.\n");
-    printout("[fatal] r_pkt: Current cell %d, target cell %d.\n", pkt.where, snext);
-    std::abort();
-  }
-
-  if (sdist > globals::max_path_step) {
-    sdist = globals::max_path_step;
-    snext = pkt.where;
-  }
 
   // At present there is no scattering/destruction process so all that needs to
   // happen is that we determine whether the packet reaches the boundary during the timestep.