From 7a7f75587c5704b1d885889b6b5408ca2d52e63a Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Fri, 7 Mar 2025 20:35:45 +0000 Subject: [PATCH] Remove last_cross --- gammapkt.cc | 11 ++--- grid.cc | 126 +++++++++++++++++++++++----------------------------- grid.h | 2 +- packet.cc | 8 +--- packet.h | 7 ++- rpkt.cc | 4 +- vpkt.cc | 4 +- 7 files changed, 65 insertions(+), 97 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 c49f50d3d..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) @@ -2432,58 +2430,41 @@ auto get_totmassradionuclide(const int z, const int a) -> double { const auto negdirections = std::array{COORD0_MIN, COORD1_MIN, COORD2_MIN}; const auto posdirections = std::array{COORD0_MAX, COORD1_MAX, COORD2_MAX}; - - 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]; - enum cell_boundary const invdirection = flip == 0 ? posdirections[d] : negdirections[d]; - const int cellindexstride = - flip != 0 ? -grid::get_coordcellindexincrement(d) : grid::get_coordcellindexincrement(d); - - bool isoutside_thisside = false; - double delta = 0.; - 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 - 10.); // 10 cm accuracy tolerance - } else { - // packet pos above max - const double boundaryposmax = cellcoordmax[d] / globals::tmin * tstart; - delta = pktposgridcoord[d] - boundaryposmax; - isoutside_thisside = pktposgridcoord[d] > (boundaryposmax + 10.); - } - - if (isoutside_thisside && (last_cross != direction)) { - // for (int d2 = 0; d2 < ndim; d2++) - const int d2 = d; - { - printout( - "[warning] packet outside coord %d %c%c boundary of cell %d. vel %g initpos %g " - "cellcoordmin %g, cellcoordmax %g\n", - d, flip != 0 ? '-' : '+', grid::coordlabel[d], cellindex, pktvelgridcoord[d2], pktposgridcoord[d2], - grid::get_cellcoordmin(cellindex, d2) / globals::tmin * tstart, - cellcoordmax[d2] / globals::tmin * tstart); + if constexpr (TESTMODE) { + 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.; + 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 - 10.); // 10 cm accuracy tolerance + } else { + // packet pos above max + const double boundaryposmax = cellcoordmax[d] / globals::tmin * tstart; + delta = pktposgridcoord[d] - boundaryposmax; + isoutside_thisside = pktposgridcoord[d] > (boundaryposmax + 10.); } - printout("globals::tmin %g tstart %g tstart/globals::tmin %g\n", globals::tmin, tstart, tstart / globals::tmin); - printout("[warning] delta %g\n", delta); - - printout("[warning] dir [%g, %g, %g]\n", dir[0], dir[1], dir[2]); - if ((pktvelgridcoord[d] - (pktposgridcoord[d] / tstart)) > 0) { - if ((grid::get_cellcoordpointnum(cellindex, d) == (grid::ncoordgrid[d] - 1) && cellindexstride > 0) || - (grid::get_cellcoordpointnum(cellindex, d) == 0 && cellindexstride < 0)) { - printout("escaping packet\n"); - return {0., -99}; + + if (isoutside_thisside) { + // for (int d2 = 0; d2 < ndim; d2++) + const int d2 = d; + { + printout( + "[warning] packet outside coord %d %c%c boundary of cell %d. vel %g initpos %g " + "cellcoordmin %g, cellcoordmax %g\n", + d, flip != 0 ? '-' : '+', grid::coordlabel[d], cellindex, pktvelgridcoord[d2], pktposgridcoord[d2], + grid::get_cellcoordmin(cellindex, d2) / globals::tmin * tstart, + cellcoordmax[d2] / globals::tmin * tstart); } - const int snext = cellindex + cellindexstride; - *pkt_last_cross = invdirection; - printout("[warning] swapping packet cellindex from %d to %d and setting last_cross to %d\n", cellindex, snext, - *pkt_last_cross); - return {0., snext}; + printout("globals::tmin %g tstart %g tstart/globals::tmin %g\n", globals::tmin, tstart, + tstart / globals::tmin); + printout("[warning] delta %g\n", delta); + + printout("[warning] dir [%g, %g, %g]\n", dir[0], dir[1], dir[2]); } - printout("pretending last_cross is %d\n", direction); - last_cross = direction; } } } @@ -2495,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; @@ -2506,14 +2486,16 @@ 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); + if ((pktvelgridcoord[1] * tstart) > pktposgridcoord[1]) { + d_coordminboundary[1] = -1; + } else { + d_coordmaxboundary[1] = -1; + } + } else if constexpr (GRID_TYPE == GridType::CARTESIAN3D) { // There are six possible boundary crossings. Each of the three // cartesian coordinates may be taken in turn. For x, the packet @@ -2528,17 +2510,22 @@ auto get_totmassradionuclide(const int z, const int a) -> double { // boundary, regardless of direction. for (int d = 0; d < 3; d++) { - const double t_coordminboundary = - ((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) / - (grid::get_cellcoordmin(cellindex, d) - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) - - tstart; - d_coordminboundary[d] = CLIGHT_PROP * t_coordminboundary; + if ((pktvelgridcoord[d] * tstart) > pktposgridcoord[d]) { + d_coordminboundary[d] = -1; + const double t_coordmaxboundary = ((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) / + (cellcoordmax[d] - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) - + tstart; - 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; + } else { + d_coordmaxboundary[d] = -1; - d_coordmaxboundary[d] = CLIGHT_PROP * t_coordmaxboundary; + const double t_coordminboundary = + ((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) / + (grid::get_cellcoordmin(cellindex, d) - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) - + tstart; + d_coordminboundary[d] = CLIGHT_PROP * t_coordminboundary; + } } } else { assert_always(false); @@ -2550,25 +2537,23 @@ 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)) { snext = -99; } else { - *pkt_last_cross = choice; snext = cellindex + grid::get_coordcellindexincrement(d); } } // 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) { snext = -99; } else { - *pkt_last_cross = choice; snext = cellindex - grid::get_coordcellindexincrement(d); } } @@ -2580,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..57525bb55 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) { @@ -166,7 +165,7 @@ void packet_init(Packet *pkt) void write_packets(const char filename[], const Packet *const pkt) { FILE *packets_file = fopen_required(filename, "w"); fprintf(packets_file, - "#number where type_id posx posy posz dirx diry dirz last_cross tdecay e_cmf e_rf nu_cmf nu_rf " + "#number where type_id posx posy posz dirx diry dirz tdecay e_cmf e_rf nu_cmf nu_rf " "escape_type_id escape_time next_trans last_event emissiontype trueemissiontype " "em_posx em_posy em_posz absorption_type absorption_freq nscatterings em_time absorptiondirx absorptiondiry " "absorptiondirz stokes1 stokes2 stokes3 pol_dirx pol_diry pol_dirz originated_from_positron " @@ -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()) {