From 98d3a7f8e64f9f833edc1842748bb04d0f081c7d Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Fri, 7 Mar 2025 21:57:21 +0000 Subject: [PATCH] Remove last_cross --- gammapkt.cc | 11 +++-------- grid.cc | 16 ++-------------- grid.h | 2 +- packet.cc | 6 ------ packet.h | 7 +++---- rpkt.cc | 4 +--- vpkt.cc | 4 +--- 7 files changed, 11 insertions(+), 39 deletions(-) diff --git a/gammapkt.cc b/gammapkt.cc index 7e8d95b32..0754fc995 100644 --- a/gammapkt.cc +++ b/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.}; diff --git a/grid.cc b/grid.cc index 5ad7f5870..0a6d52ddd 100644 --- a/grid.cc +++ b/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 &dir, const std::array &pos, const double tstart, - const int cellindex, enum cell_boundary *pkt_last_cross) - -> std::tuple { + const int cellindex) -> std::tuple { 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{-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]); } diff --git a/grid.h b/grid.h index e02c8f5e3..15d448d5c 100644 --- a/grid.h +++ b/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 &dir, const std::array &pos, double tstart, - int cellindex, enum cell_boundary *pkt_last_cross) -> std::tuple; + int cellindex) -> std::tuple; void calculate_kappagrey(); diff --git a/packet.cc b/packet.cc index e42701f33..0317302d7 100644 --- a/packet.cc +++ b/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(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(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(last_cross_in); - ssline >> pkt[i].tdecay; ssline >> pkt[i].e_cmf >> pkt[i].e_rf >> pkt[i].nu_cmf >> pkt[i].nu_rf; diff --git a/packet.h b/packet.h index 6962558a1..81897f022 100644 --- a/packet.h +++ b/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 pos{}; // Position of the packet (x,y,z). diff --git a/rpkt.cc b/rpkt.cc index 19506ca08..c8f523064 100644 --- a/rpkt.cc +++ b/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 diff --git a/vpkt.cc b/vpkt.cc index 938559b66..568338a6e 100644 --- a/vpkt.cc +++ b/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()) {