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 8, 2025
1 parent 307060d commit b48c2b8
Show file tree
Hide file tree
Showing 7 changed files with 75 additions and 106 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
145 changes: 65 additions & 80 deletions grid.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1317,9 +1317,12 @@ auto get_coordboundary_distances_cylindrical2d(const std::array<double, 3> &pkt_
d_coordminboundary[0] = -1;
// don't try to calculate the intersection if the inner radius is zero
if (r_inner > 0) {
// calculate the distance in the xy plane to the inner boundary
const double d_rcyl_coordminboundary = expanding_shell_intersection(posnoz, dirnoz, xyspeed, r_inner, true, tstart);
if (d_rcyl_coordminboundary >= 0) {
// how far did the packet travel in the z direction during this time?
const double d_z_coordminboundary = d_rcyl_coordminboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP;
// get distance from two perpendicular components
d_coordminboundary[0] = std::sqrt((d_rcyl_coordminboundary * d_rcyl_coordminboundary) +
(d_z_coordminboundary * d_z_coordminboundary));
}
Expand All @@ -1334,17 +1337,24 @@ auto get_coordboundary_distances_cylindrical2d(const std::array<double, 3> &pkt_
std::sqrt((d_rcyl_coordmaxboundary * d_rcyl_coordmaxboundary) + (d_z_coordmaxboundary * d_z_coordmaxboundary));
}

// z boundaries are the same as Cartesian
const double t_zcoordminboundary =
((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) /
((grid::get_cellcoordmin(cellindex, 1)) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) -
tstart;
d_coordminboundary[1] = CLIGHT_PROP * t_zcoordminboundary;
// 2D cylindrical Z boundaries are the same as Cartesian

const double t_zcoordmaxboundary = ((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) /
((cellcoordmax[1]) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) -
tstart;
d_coordmaxboundary[1] = CLIGHT_PROP * t_zcoordmaxboundary;
if ((pktvelgridcoord[1] * tstart) > pktposgridcoord[1]) {
d_coordminboundary[1] = -1;

const double t_zcoordmaxboundary = ((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) /
((cellcoordmax[1]) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) -
tstart;
d_coordmaxboundary[1] = CLIGHT_PROP * t_zcoordmaxboundary;
} else {
const double t_zcoordminboundary =
((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) /
((grid::get_cellcoordmin(cellindex, 1)) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) -
tstart;
d_coordminboundary[1] = CLIGHT_PROP * t_zcoordminboundary;

d_coordmaxboundary[1] = -1;
}

return {d_coordminboundary, d_coordmaxboundary};
}
Expand Down Expand Up @@ -2390,15 +2400,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 +2440,38 @@ 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};
if constexpr (TESTMODE || true) {
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.;
const bool pos_component_vel_relative_to_flow = (pktvelgridcoord[d] * tstart) > pktposgridcoord[d];
if (flip != 0) {
// packet pos below min
const double boundaryposmin = grid::get_cellcoordmin(cellindex, d) / globals::tmin * tstart;
delta = pktposgridcoord[d] - boundaryposmin;
isoutside_thisside = !pos_component_vel_relative_to_flow && (pktposgridcoord[d] < (boundaryposmin - 10.));
} else {
// packet pos above max
const double boundaryposmax = cellcoordmax[d] / globals::tmin * tstart;
delta = pktposgridcoord[d] - boundaryposmax;
isoutside_thisside = pos_component_vel_relative_to_flow && (pktposgridcoord[d] > (boundaryposmax + 10.));
}

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;
{
if (isoutside_thisside) {
printout(
"[warning] packet outside coord %d %c%c boundary of cell %d. vel %g initpos %g "
"[ERROR] 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);
}
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};
}
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};
d, flip != 0 ? '-' : '+', grid::coordlabel[d], cellindex, pktvelgridcoord[d], pktposgridcoord[d],
grid::get_cellcoordmin(cellindex, d) / globals::tmin * tstart, cellcoordmax[d] / globals::tmin * tstart);
printout("globals::tmin %g tstart %g tstart/globals::tmin %g\n", globals::tmin, tstart,
tstart / globals::tmin);
printout(" delta %g\n", delta);

printout("packet dir [%g, %g, %g]\n", dir[0], dir[1], dir[2]);
assert_always(false);
}
printout("pretending last_cross is %d\n", direction);
last_cross = direction;
}
}
}
Expand All @@ -2495,7 +2483,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,10 +2493,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);
Expand All @@ -2528,17 +2511,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 +2538,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 +2566,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
8 changes: 1 addition & 7 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 @@ -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 "
Expand All @@ -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
Loading

0 comments on commit b48c2b8

Please sign in to comment.