From 477b25b7775376a89ac315b9ffac98b5b2724b27 Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Fri, 7 Mar 2025 22:43:22 +0000 Subject: [PATCH] Partial revert --- gammapkt.cc | 11 ++++++++--- grid.cc | 16 ++++++++++++++-- grid.h | 2 +- packet.cc | 6 ++++++ packet.h | 7 ++++--- rpkt.cc | 4 +++- vpkt.cc | 4 +++- 7 files changed, 39 insertions(+), 11 deletions(-) diff --git a/gammapkt.cc b/gammapkt.cc index 0754fc995..7e8d95b32 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -443,6 +443,8 @@ 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) { @@ -711,6 +713,7 @@ void pair_prod(Packet &pkt) { pkt.e_rf = pkt.e_cmf / dopplerfactor; pkt.type = TYPE_GAMMA; + pkt.last_cross = BOUNDARY_NONE; } } @@ -726,7 +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. - const auto [sdist, snext] = grid::boundary_distance(pkt.dir, pkt.pos, pkt.prop_time, 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. @@ -863,7 +866,8 @@ 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); + 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 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()) { @@ -925,7 +929,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); + grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, &pkt_copy.last_cross); 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()) { @@ -1017,6 +1021,7 @@ __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 0a6d52ddd..5ad7f5870 100644 --- a/grid.cc +++ b/grid.cc @@ -2390,13 +2390,15 @@ 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) -> std::tuple { + const int cellindex, enum cell_boundary *pkt_last_cross) + -> 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) @@ -2434,6 +2436,8 @@ 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) { @@ -2448,7 +2452,7 @@ auto get_totmassradionuclide(const int z, const int a) -> double { isoutside_thisside = pktposgridcoord[d] > (boundaryposmax + 10.); } - if (isoutside_thisside) { + if (isoutside_thisside && (last_cross != direction)) { // for (int d2 = 0; d2 < ndim; d2++) const int d2 = d; { @@ -2476,6 +2480,7 @@ 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; @@ -2486,6 +2491,10 @@ 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); @@ -2543,6 +2552,7 @@ 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); } } @@ -2554,6 +2564,7 @@ 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); } } @@ -2565,6 +2576,7 @@ 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 15d448d5c..e02c8f5e3 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) -> std::tuple; + int cellindex, enum cell_boundary *pkt_last_cross) -> std::tuple; void calculate_kappagrey(); diff --git a/packet.cc b/packet.cc index 0317302d7..e42701f33 100644 --- a/packet.cc +++ b/packet.cc @@ -32,6 +32,7 @@ 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) { @@ -176,6 +177,7 @@ 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); @@ -277,6 +279,10 @@ 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 81897f022..6962558a1 100644 --- a/packet.h +++ b/packet.h @@ -40,9 +40,10 @@ 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 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. 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 c8f523064..19506ca08 100644 --- a/rpkt.cc +++ b/rpkt.cc @@ -292,6 +292,7 @@ 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); @@ -625,7 +626,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); + 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); @@ -956,6 +957,7 @@ __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 568338a6e..938559b66 100644 --- a/vpkt.cc +++ b/vpkt.cc @@ -168,6 +168,7 @@ 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; @@ -253,7 +254,8 @@ 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); + const auto [sdist, snext] = + grid::boundary_distance(vpkt.dir, vpkt.pos, vpkt.prop_time, vpkt.where, &vpkt.last_cross); const double s_cont = sdist * t_current * t_current * t_current / (t_future * t_future * t_future); if (mgi == grid::get_npts_model()) {