Skip to content

Commit

Permalink
Remove last_cross
Browse files Browse the repository at this point in the history
  • Loading branch information
@lukeshingles
lukeshingles committed Mar 7, 2025
1 parent 95ec30c commit d9c9544
Show file tree
Hide file tree
Showing 7 changed files with 64 additions and 96 deletions.
11 changes: 3 additions & 8 deletions gammapkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
}
}

Expand All @@ -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.
Expand Down Expand Up @@ -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()) {
Expand Down Expand Up @@ -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()) {
Expand Down Expand Up @@ -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.};
Expand Down
126 changes: 55 additions & 71 deletions grid.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -2432,58 +2430,41 @@ auto get_totmassradionuclide(const int z, const int a) -> double {

const auto negdirections = std::array<enum cell_boundary, 3>{COORD0_MIN, COORD1_MIN, COORD2_MIN};
const auto posdirections = std::array<enum cell_boundary, 3>{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;
}
}
}
Expand All @@ -2495,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;
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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);
}
}
Expand All @@ -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]);
}
Expand Down
2 changes: 1 addition & 1 deletion grid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down
6 changes: 0 additions & 6 deletions packet.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down
7 changes: 3 additions & 4 deletions packet.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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).
Expand Down
4 changes: 1 addition & 3 deletions rpkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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

Expand Down
4 changes: 1 addition & 3 deletions vpkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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()) {
Expand Down

0 comments on commit d9c9544

Please sign in to comment.