Remove last_cross

gammapkt.cc

@@ -443,8 +443,6 @@ void compton_scatter(Packet &pkt) {
     const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time);
     pkt.nu_rf = pkt.nu_cmf / dopplerfactor;
     pkt.e_rf = pkt.e_cmf / dopplerfactor;
-
-    pkt.last_cross = BOUNDARY_NONE;  // allow it to re-cross a boundary
   } else {
     // energy loss of the gamma becomes energy of the electron (needed to calculate time-dependent thermalisation rate)
     if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
@@ -713,7 +711,6 @@ void pair_prod(Packet &pkt) {
     pkt.e_rf = pkt.e_cmf / dopplerfactor;
 
     pkt.type = TYPE_GAMMA;
-    pkt.last_cross = BOUNDARY_NONE;
   }
 }
 
@@ -729,7 +726,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.
 
-  const 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);
 
   // Now consider the scattering/destruction processes.
   // Compton scattering - need to determine the scattering co-efficient.
@@ -866,8 +863,7 @@ void wollaeger_thermalisation(Packet &pkt) {
   bool end_packet = false;
   while (!end_packet) {
     // distance to the next cell
-    const auto [sdist, snext] =
-        grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, &pkt_copy.last_cross);
+    const auto [sdist, snext] = grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where);
     const double s_cont = sdist * t_current * t_current * t_current / std::pow(pkt_copy.prop_time, 3);
     const int mgi = grid::get_propcell_modelgridindex(pkt_copy.where);
     if (mgi != grid::get_npts_model()) {
@@ -929,7 +925,7 @@ void guttman_thermalisation(Packet &pkt) {
     while (!end_packet) {
       // distance to the next cell
       const auto [sdist, snext] =
-          grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, &pkt_copy.last_cross);
+          grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where);
       const double s_cont = sdist * std::pow(t, 3.) / std::pow(pkt_copy.prop_time, 3.);
       const int mgi = grid::get_propcell_modelgridindex(pkt_copy.where);
       if (mgi != grid::get_npts_model()) {
@@ -1021,7 +1017,6 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) {
   pkt.e_rf = pkt.e_cmf / dopplerfactor;
 
   pkt.type = TYPE_GAMMA;
-  pkt.last_cross = BOUNDARY_NONE;
 
   // initialise polarisation information
   pkt.stokes = {1., 0., 0.};

grid.cc

@@ -2390,15 +2390,13 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
 // compute distance to a cell boundary.
 [[nodiscard]] __host__ __device__ auto boundary_distance(const std::array<double, 3> &dir,
                                                          const std::array<double, 3> &pos, const double tstart,
-                                                         const int cellindex, enum cell_boundary *pkt_last_cross)
-    -> std::tuple<double, int> {
+                                                         const int cellindex) -> std::tuple<double, int> {
   if constexpr (FORCE_SPHERICAL_ESCAPE_SURFACE) {
     if (get_cell_r_inner(cellindex) > globals::vmax * globals::tmin) {
       return {0., -99};
     }
   }
 
-  auto last_cross = *pkt_last_cross;
   // d is used to loop over the coordinate indicies 0,1,2 for x,y,z
 
   // the following vector are in grid coordinates, so either x,y,z (3D) or r (1D), or r_xy, z (2D)
@@ -2436,8 +2434,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
     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++) {
-        enum cell_boundary const direction = flip != 0 ? posdirections[d] : negdirections[d];
-
         bool isoutside_thisside = false;
         double delta = 0.;
         if (flip != 0) {
@@ -2452,7 +2448,7 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
           isoutside_thisside = pktposgridcoord[d] > (boundaryposmax + 10.);
         }
 
-        if (isoutside_thisside && (last_cross != direction)) {
+        if (isoutside_thisside) {
           // for (int d2 = 0; d2 < ndim; d2++)
           const int d2 = d;
           {
@@ -2480,7 +2476,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
   auto d_coordminboundary = std::array<double, 3>{-1};
 
   if constexpr (GRID_TYPE == GridType::SPHERICAL1D) {
-    last_cross = BOUNDARY_NONE;  // handle this separately by setting d_inner and d_outer negative for invalid direction
     const double speed = vec_len(dir) * CLIGHT_PROP;  // just in case dir is not normalised
 
     const double r_inner = grid::get_cellcoordmin(cellindex, 0) * tstart / globals::tmin;
@@ -2491,10 +2486,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
 
   } else if constexpr (GRID_TYPE == GridType::CYLINDRICAL2D) {
     // coordinate 0 is radius in x-y plane, coord 1 is z
-    if (last_cross == COORD0_MIN || last_cross == COORD0_MAX) {
-      last_cross =
-          BOUNDARY_NONE;  // handle this separately by setting d_inner and d_outer negative for invalid direction
-    }
 
     std::tie(d_coordminboundary, d_coordmaxboundary) = get_coordboundary_distances_cylindrical2d(
         pos, dir, pktposgridcoord, pktvelgridcoord, cellindex, tstart, cellcoordmax);
@@ -2552,7 +2543,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
       if (grid::get_cellcoordpointnum(cellindex, d) == (grid::ncoordgrid[d] - 1)) {
         snext = -99;
       } else {
-        *pkt_last_cross = choice;
         snext = cellindex + grid::get_coordcellindexincrement(d);
       }
     }
@@ -2564,7 +2554,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
       if (grid::get_cellcoordpointnum(cellindex, d) == 0) {
         snext = -99;
       } else {
-        *pkt_last_cross = choice;
         snext = cellindex - grid::get_coordcellindexincrement(d);
       }
     }
@@ -2576,7 +2565,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
       printout("packet cell %d\n", cellindex);
       printout("choice %d\n", choice);
       printout("globals::tmin %g tstart %g\n", globals::tmin, tstart);
-      printout("last_cross %d\n", last_cross);
       for (int d2 = 0; d2 < 3; d2++) {
         printout("coord %d: initpos %g dir %g\n", d2, pos[d2], dir[d2]);
       }

grid.h

@@ -106,7 +106,7 @@ void write_grid_restart_data(int timestep);
 [[nodiscard]] auto get_ndo_nonempty(int rank) -> int;
 [[nodiscard]] auto get_totmassradionuclide(int z, int a) -> double;
 [[nodiscard]] auto boundary_distance(const std::array<double, 3> &dir, const std::array<double, 3> &pos, double tstart,
-                                     int cellindex, enum cell_boundary *pkt_last_cross) -> std::tuple<double, int>;
+                                     int cellindex) -> std::tuple<double, int>;
 
 void calculate_kappagrey();
 

packet.cc

@@ -32,7 +32,6 @@ void place_pellet(const double e0, const int cellindex, const int pktnumber, Pac
   pkt.where = cellindex;
   pkt.number = pktnumber;  // record the packets number for debugging
   pkt.prop_time = globals::tmin;
-  // pkt.last_cross = BOUNDARY_NONE;
   pkt.originated_from_particlenotgamma = false;
 
   if constexpr (GRID_TYPE == GridType::SPHERICAL1D) {
@@ -177,7 +176,6 @@ void write_packets(const char filename[], const Packet *const pkt) {
     fprintf(packets_file, "%d ", static_cast<int>(pkt[i].type));
     fprintf(packets_file, "%lg %lg %lg ", pkt[i].pos[0], pkt[i].pos[1], pkt[i].pos[2]);
     fprintf(packets_file, "%lg %lg %lg ", pkt[i].dir[0], pkt[i].dir[1], pkt[i].dir[2]);
-    fprintf(packets_file, "%d ", static_cast<int>(pkt[i].last_cross));
     fprintf(packets_file, "%g ", pkt[i].tdecay);
     fprintf(packets_file, "%g ", pkt[i].e_cmf);
     fprintf(packets_file, "%g ", pkt[i].e_rf);
@@ -279,10 +277,6 @@ void read_packets(const char filename[], Packet *pkt) {
 
     ssline >> pkt[i].dir[0] >> pkt[i].dir[1] >> pkt[i].dir[2];
 
-    int last_cross_in = 0;
-    ssline >> last_cross_in;
-    pkt[i].last_cross = static_cast<enum cell_boundary>(last_cross_in);
-
     ssline >> pkt[i].tdecay;
 
     ssline >> pkt[i].e_cmf >> pkt[i].e_rf >> pkt[i].nu_cmf >> pkt[i].nu_rf;

packet.h

@@ -40,10 +40,9 @@ enum cell_boundary : int {
 };
 
 struct Packet {
-  enum packet_type type {};  // type of packet (k-, r-, etc.)
-  double prop_time{-1.};     // internal clock to track how far in time the packet has been propagated
-  int where{-1};             // The propagation grid cell that the packet is in.
-  enum cell_boundary last_cross { BOUNDARY_NONE };  // To avoid rounding errors on cell crossing.
+  enum packet_type type {};           // type of packet (k-, r-, etc.)
+  double prop_time{-1.};              // internal clock to track how far in time the packet has been propagated
+  int where{-1};                      // The propagation grid cell that the packet is in.
   int nscatterings{0};                // records number of electron scatterings a r-pkt undergone since it was emitted
   int last_event{0};                  // debug: stores information about the packets history
   std::array<double, 3> pos{};        // Position of the packet (x,y,z).

rpkt.cc

@@ -292,7 +292,6 @@ auto get_event_expansion_opacity(
 void electron_scatter_rpkt(Packet &pkt) {
   // now make the packet a r-pkt and set further flags
   pkt.type = TYPE_RPKT;
-  pkt.last_cross = BOUNDARY_NONE;  // allow all further cell crossings
 
   const auto vel_vec = get_velocity(pkt.pos, pkt.prop_time);
 
@@ -626,7 +625,7 @@ 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.
-  const 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);
 
   if (sdist == 0) {
     grid::change_cell(pkt, snext);
@@ -957,7 +956,6 @@ __host__ __device__ void do_rpkt(Packet &pkt, const double t2) {
 // make the packet an r-pkt and set further flags
 __host__ __device__ void emit_rpkt(Packet &pkt) {
   pkt.type = TYPE_RPKT;
-  pkt.last_cross = BOUNDARY_NONE;  // allow all further cell crossings
 
   // Need to assign a new direction. Assume isotropic emission in the cmf
 

vpkt.cc

@@ -168,7 +168,6 @@ auto rlc_emiss_vpkt(const Packet &pkt, const double t_current, const double t_ar
   vpkt.nu_rf = nu_rf;
   vpkt.e_rf = e_rf;
   vpkt.dir = obsdir;
-  vpkt.last_cross = BOUNDARY_NONE;
 
   bool end_packet = false;
   double t_future = t_current;
@@ -254,8 +253,7 @@ auto rlc_emiss_vpkt(const Packet &pkt, const double t_current, const double t_ar
 
   while (!end_packet) {
     // distance to the next cell
-    const auto [sdist, snext] =
-        grid::boundary_distance(vpkt.dir, vpkt.pos, vpkt.prop_time, vpkt.where, &vpkt.last_cross);
+    const auto [sdist, snext] = grid::boundary_distance(vpkt.dir, vpkt.pos, vpkt.prop_time, vpkt.where);
     const double s_cont = sdist * t_current * t_current * t_current / (t_future * t_future * t_future);
 
     if (mgi == grid::get_npts_model()) {