Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Include neglect expansion option in grid boundary dist calc #101

Open
wants to merge 2 commits into
base: develop
Choose a base branch
from
+34 −26
Open

Include neglect expansion option in grid boundary dist calc #101

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d482f5e
Include neglect expansion option in grid boundary dist calc
gleck97 Aug 7, 2024
bc02401
move_pkt_withtime: Use template parameter for CLIGHT_ASSUMED
lukeshingles Aug 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 7 additions & 6 deletions gammapkt.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1000,15 +1000,15 @@ 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,
&pkt_copy.last_cross, true);
const double s_cont = sdist * t_current * t_current * t_current / std::pow(pkt_copy.prop_time, 3);
const int mgi = grid::get_cell_modelgridindex(pkt_copy.where);
if (mgi != grid::get_npts_model()) {
tau += grid::get_rho(mgi) * s_cont * mean_gamma_opac; // contribution to the integral
}
// move packet copy now
move_pkt_withtime(pkt_copy, sdist);
move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist);

grid::change_cell(pkt_copy, snext);
end_packet = (pkt_copy.type == TYPE_ESCAPE);
Expand Down Expand Up @@ -1061,15 +1061,16 @@ void guttman_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, &pkt_copy.last_cross, true);
const double s_cont = sdist * std::pow(t, 3.) / std::pow(pkt_copy.prop_time, 3.);
std::cout << "t: " << t << " s_cont: " << s_cont << "\n";
const int mgi = grid::get_cell_modelgridindex(pkt_copy.where);
if (mgi != grid::get_npts_model()) {
column_densities[i] += grid::get_rho_tmin(mgi) * s_cont; // contribution to the integral
}
// move packet copy now
move_pkt_withtime(pkt_copy, sdist);
move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist);

grid::change_cell(pkt_copy, snext);
end_packet = (pkt_copy.type == TYPE_ESCAPE);
Expand Down
37 changes: 20 additions & 17 deletions grid.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2465,18 +2465,19 @@ template <size_t S1>
static auto get_coordboundary_distances_cylindrical2d(
const std::array<double, 3> pkt_pos, const std::array<double, 3> pkt_dir,
const std::array<double, 3> pktposgridcoord, const std::array<double, 3> pktvelgridcoord, int cellindex,
const double tstart,
const std::array<double, 3> cellcoordmax) -> std::tuple<std::array<double, 3>, std::array<double, 3>> {
const double tstart, const std::array<double, 3> cellcoordmax,
bool neglect_exp) -> std::tuple<std::array<double, 3>, std::array<double, 3>> {
// to get the cylindrical intersection, get the spherical intersection with Z components set to zero, and the
// propagation speed set to the xy component of the 3-velocity
double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP;

std::array<double, 3> d_coordminboundary{};
std::array<double, 3> d_coordmaxboundary{};

const std::array<double, 2> posnoz = {pkt_pos[0], pkt_pos[1]};

const double dirxylen = std::sqrt((pkt_dir[0] * pkt_dir[0]) + (pkt_dir[1] * pkt_dir[1]));
const double xyspeed = dirxylen * CLIGHT_PROP; // r_cyl component of velocity
const double xyspeed = dirxylen * CLIGHT_PROP_DIST_CALC; // r_cyl component of velocity

// make a normalised direction vector in the xy plane
const std::array<double, 2> dirnoz = {pkt_dir[0] / dirxylen, pkt_dir[1] / dirxylen};
Expand All @@ -2487,7 +2488,7 @@ static auto get_coordboundary_distances_cylindrical2d(
if (r_inner > 0) {
const double d_rcyl_coordminboundary = expanding_shell_intersection(posnoz, dirnoz, xyspeed, r_inner, true, tstart);
if (d_rcyl_coordminboundary >= 0) {
const double d_z_coordminboundary = d_rcyl_coordminboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP;
const double d_z_coordminboundary = d_rcyl_coordminboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP_DIST_CALC;
d_coordminboundary[0] =
std::sqrt(d_rcyl_coordminboundary * d_rcyl_coordminboundary + d_z_coordminboundary * d_z_coordminboundary);
}
Expand All @@ -2497,7 +2498,7 @@ static auto get_coordboundary_distances_cylindrical2d(
const double d_rcyl_coordmaxboundary = expanding_shell_intersection(posnoz, dirnoz, xyspeed, r_outer, false, tstart);
d_coordmaxboundary[0] = -1;
if (d_rcyl_coordmaxboundary >= 0) {
const double d_z_coordmaxboundary = d_rcyl_coordmaxboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP;
const double d_z_coordmaxboundary = d_rcyl_coordmaxboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP_DIST_CALC;
d_coordmaxboundary[0] =
std::sqrt(d_rcyl_coordmaxboundary * d_rcyl_coordmaxboundary + d_z_coordmaxboundary * d_z_coordmaxboundary);
}
Expand All @@ -2507,21 +2508,23 @@ static auto get_coordboundary_distances_cylindrical2d(
((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_coordminboundary[1] = CLIGHT_PROP_DIST_CALC * t_zcoordminboundary;

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;
d_coordmaxboundary[1] = CLIGHT_PROP_DIST_CALC * t_zcoordmaxboundary;

return {d_coordminboundary, d_coordmaxboundary};
}

[[nodiscard]] 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 double tstart, const int cellindex, enum cell_boundary *pkt_last_cross,
bool neglect_exp) -> std::tuple<double, int>
/// Basic routine to compute distance to a cell boundary.
{
double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP;

if constexpr (FORCE_SPHERICAL_ESCAPE_SURFACE) {
if (get_cell_r_inner(cellindex) > globals::vmax * globals::tmin) {
return {0., -99};
Expand All @@ -2545,18 +2548,18 @@ static auto get_coordboundary_distances_cylindrical2d(
if constexpr (GRID_TYPE == GRID_CARTESIAN3D) {
// keep xyz Cartesian coordinates
for (int d = 0; d < ndim; d++) {
pktvelgridcoord[d] = dir[d] * CLIGHT_PROP;
pktvelgridcoord[d] = dir[d] * CLIGHT_PROP_DIST_CALC;
}
} else if constexpr (GRID_TYPE == GRID_CYLINDRICAL2D) {
// xy plane radial velocity
pktvelgridcoord[0] = (pos[0] * dir[0] + pos[1] * dir[1]) / pktposgridcoord[0] * CLIGHT_PROP;
pktvelgridcoord[0] = (pos[0] * dir[0] + pos[1] * dir[1]) / pktposgridcoord[0] * CLIGHT_PROP_DIST_CALC;

// second cylindrical coordinate is z
pktvelgridcoord[1] = dir[2] * CLIGHT_PROP;
pktvelgridcoord[1] = dir[2] * CLIGHT_PROP_DIST_CALC;

} else if constexpr (GRID_TYPE == GRID_SPHERICAL1D) {
// the only coordinate is radius from the origin
pktvelgridcoord[0] = dot(pos, dir) / pktposgridcoord[0] * CLIGHT_PROP; // radial velocity
pktvelgridcoord[0] = dot(pos, dir) / pktposgridcoord[0] * CLIGHT_PROP_DIST_CALC; // radial velocity
} else {
assert_always(false);
}
Expand Down Expand Up @@ -2634,7 +2637,7 @@ static auto get_coordboundary_distances_cylindrical2d(
std::array<double, 3>{-1}; // distance to reach the cell's lower boundary on each coordinate
if constexpr (GRID_TYPE == GRID_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 speed = vec_len(dir) * CLIGHT_PROP_DIST_CALC; // just in case dir is not normalised

const double r_inner = grid::get_cellcoordmin(cellindex, 0) * tstart / globals::tmin;
d_coordminboundary[0] = (r_inner > 0.) ? expanding_shell_intersection(pos, dir, speed, r_inner, true, tstart) : -1.;
Expand All @@ -2650,7 +2653,7 @@ static auto get_coordboundary_distances_cylindrical2d(
}

std::tie(d_coordminboundary, d_coordmaxboundary) = get_coordboundary_distances_cylindrical2d(
pos, dir, pktposgridcoord, pktvelgridcoord, cellindex, tstart, cellcoordmax);
pos, dir, pktposgridcoord, pktvelgridcoord, cellindex, tstart, cellcoordmax, neglect_exp);

} else if constexpr (GRID_TYPE == GRID_CARTESIAN3D) {
// There are six possible boundary crossings. Each of the three
Expand All @@ -2670,13 +2673,13 @@ static auto get_coordboundary_distances_cylindrical2d(
((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) /
(grid::get_cellcoordmin(cellindex, d) - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) -
tstart;
d_coordminboundary[d] = CLIGHT_PROP * t_coordminboundary;
d_coordminboundary[d] = CLIGHT_PROP_DIST_CALC * t_coordminboundary;

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;
d_coordmaxboundary[d] = CLIGHT_PROP_DIST_CALC * t_coordmaxboundary;
}
} else {
assert_always(false);
Expand Down
3 changes: 2 additions & 1 deletion grid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,8 @@ 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(std::array<double, 3> dir, std::array<double, 3> pos, double tstart, int cellindex,
enum cell_boundary *pkt_last_cross) -> std::tuple<double, int>;
enum cell_boundary *pkt_last_cross,
bool neglect_exp = false) -> std::tuple<double, int>;

[[nodiscard]] inline auto get_elem_abundance(int modelgridindex, int element) -> float
// mass fraction of an element (all isotopes combined)
Expand Down
7 changes: 5 additions & 2 deletions vectors.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,7 @@ template <size_t S1, size_t S2>
return doppler_nucmf_on_nurf(dir_rf, get_velocity(pos_rf, prop_time));
}

template <double CLIGHT_ASSUMED = CLIGHT_PROP>
constexpr auto move_pkt_withtime(std::span<double, 3> pos_rf, const std::array<double, 3> dir_rf, double &prop_time,
const double nu_rf, double &nu_cmf, const double e_rf, double &e_cmf,
const double distance) -> double
Expand All @@ -150,7 +151,7 @@ constexpr auto move_pkt_withtime(std::span<double, 3> pos_rf, const std::array<d
assert_always(distance >= 0);

const double nu_cmf_old = nu_cmf;
prop_time += distance / CLIGHT_PROP;
prop_time += distance / CLIGHT_ASSUMED;

pos_rf[0] += (dir_rf[0] * distance);
pos_rf[1] += (dir_rf[1] * distance);
Expand All @@ -170,8 +171,10 @@ constexpr auto move_pkt_withtime(std::span<double, 3> pos_rf, const std::array<d
return dopplerfactor;
}

template <double CLIGHT_ASSUMED = CLIGHT_PROP>
constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
return move_pkt_withtime(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, pkt.nu_cmf, pkt.e_rf, pkt.e_cmf, distance);
return move_pkt_withtime<CLIGHT_ASSUMED>(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, pkt.nu_cmf, pkt.e_rf, pkt.e_cmf,
distance);
}

[[nodiscard]] [[gnu::const]] constexpr auto get_arrive_time(const Packet &pkt) -> double
Expand Down