From d482f5ed4cbb9828318c326fab625d7bd48e9f4a Mon Sep 17 00:00:00 2001 From: gleck97 Date: Wed, 7 Aug 2024 16:24:53 +0200 Subject: [PATCH 1/2] Include neglect expansion option in grid boundary dist calc Co-Authored-By: gleck97 <86471143+gleck97@users.noreply.github.com> --- gammapkt.cc | 15 ++++++++------- grid.cc | 45 +++++++++++++++++++++++++-------------------- grid.h | 3 ++- vectors.h | 26 ++++++++++++++------------ 4 files changed, 49 insertions(+), 40 deletions(-) diff --git a/gammapkt.cc b/gammapkt.cc index 3f49c1726..1134e30c6 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -409,8 +409,8 @@ static auto thomson_angle() -> double { return mu; } -[[nodiscard]] static auto scatter_dir(const std::array dir_in, - const double cos_theta) -> std::array +[[nodiscard]] static auto scatter_dir(const std::array dir_in, const double cos_theta) + -> std::array // Routine for scattering a direction through angle theta. { // begin with setting the direction in coordinates where original direction @@ -1001,14 +1001,14 @@ void wollaeger_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, &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(pkt_copy, sdist, true); grid::change_cell(pkt_copy, snext); end_packet = (pkt_copy.type == TYPE_ESCAPE); @@ -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(pkt_copy, sdist, true); grid::change_cell(pkt_copy, snext); end_packet = (pkt_copy.type == TYPE_ESCAPE); diff --git a/grid.cc b/grid.cc index 3ddc27b4c..b0107847c 100644 --- a/grid.cc +++ b/grid.cc @@ -2462,13 +2462,16 @@ template return -1.; } -static auto get_coordboundary_distances_cylindrical2d( - const std::array pkt_pos, const std::array pkt_dir, - const std::array pktposgridcoord, const std::array pktvelgridcoord, int cellindex, - const double tstart, - const std::array cellcoordmax) -> std::tuple, std::array> { +static auto get_coordboundary_distances_cylindrical2d(const std::array pkt_pos, + const std::array pkt_dir, + const std::array pktposgridcoord, + const std::array pktvelgridcoord, int cellindex, + const double tstart, const std::array cellcoordmax, + bool neglect_exp) + -> std::tuple, std::array> { // 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 d_coordminboundary{}; std::array d_coordmaxboundary{}; @@ -2476,7 +2479,7 @@ static auto get_coordboundary_distances_cylindrical2d( const std::array 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 dirnoz = {pkt_dir[0] / dirxylen, pkt_dir[1] / dirxylen}; @@ -2487,7 +2490,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); } @@ -2497,7 +2500,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); } @@ -2507,21 +2510,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 dir, const std::array pos, - const double tstart, const int cellindex, - enum cell_boundary *pkt_last_cross) -> std::tuple + const double tstart, const int cellindex, enum cell_boundary *pkt_last_cross, + bool neglect_exp) -> std::tuple /// 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}; @@ -2545,18 +2550,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); } @@ -2634,7 +2639,7 @@ static auto get_coordboundary_distances_cylindrical2d( std::array{-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.; @@ -2650,7 +2655,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 @@ -2670,13 +2675,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); diff --git a/grid.h b/grid.h index 72bba2793..95db53000 100644 --- a/grid.h +++ b/grid.h @@ -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 dir, std::array pos, double tstart, int cellindex, - enum cell_boundary *pkt_last_cross) -> std::tuple; + enum cell_boundary *pkt_last_cross, bool neglect_exp = false) + -> std::tuple; [[nodiscard]] inline auto get_elem_abundance(int modelgridindex, int element) -> float // mass fraction of an element (all isotopes combined) diff --git a/vectors.h b/vectors.h index df854cd74..e70be76f1 100644 --- a/vectors.h +++ b/vectors.h @@ -34,15 +34,15 @@ template } template -[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array x, - const std::array y) -> double +[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array x, const std::array y) + -> double // vector dot product { return std::inner_product(x.begin(), x.end(), y.begin(), 0.); } -[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span x, - const double t) -> std::array +[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span x, const double t) + -> std::array // Routine for getting velocity vector of the flow at a position with homologous expansion. { return std::array{x[0] / t, x[1] / t, x[2] / t}; @@ -59,8 +59,8 @@ template return std::array{vec[0] * scalefactor, vec[1] * scalefactor, vec[2] * scalefactor}; } -[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array dir1, - const std::array vel) -> std::array +[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array dir1, const std::array vel) + -> std::array // aberation of angles in special relativity // dir1: direction unit vector in frame1 // vel: velocity of frame2 relative to frame1 @@ -143,14 +143,15 @@ template constexpr auto move_pkt_withtime(std::span pos_rf, const std::array dir_rf, double &prop_time, const double nu_rf, double &nu_cmf, const double e_rf, double &e_cmf, - const double distance) -> double + const double distance, bool neglect_exp = false) -> double /// Subroutine to move a packet along a straight line (specified by current /// dir vector). The distance moved is in the rest frame. { + double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP; assert_always(distance >= 0); const double nu_cmf_old = nu_cmf; - prop_time += distance / CLIGHT_PROP; + prop_time += distance / CLIGHT_PROP_DIST_CALC; pos_rf[0] += (dir_rf[0] * distance); pos_rf[1] += (dir_rf[1] * distance); @@ -170,8 +171,9 @@ constexpr auto move_pkt_withtime(std::span pos_rf, const std::array 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); +constexpr auto move_pkt_withtime(Packet &pkt, const double distance, bool neglect_exp = false) -> 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, + neglect_exp); } [[nodiscard]] [[gnu::const]] constexpr auto get_arrive_time(const Packet &pkt) -> double @@ -308,8 +310,8 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double { } // Routine to transform the Stokes Parameters from RF to CMF -constexpr auto frame_transform(const std::array n_rf, double *Q, double *U, - const std::array v) -> std::array { +constexpr auto frame_transform(const std::array n_rf, double *Q, double *U, const std::array v) + -> std::array { // Meridian frame in the RF const auto [ref1_rf, ref2_rf] = meridian(n_rf); From bc024016cedaddfdf45cb49166a5626c88f054d3 Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Fri, 16 Aug 2024 12:00:27 +0100 Subject: [PATCH 2/2] move_pkt_withtime: Use template parameter for CLIGHT_ASSUMED --- gammapkt.cc | 12 ++++++------ grid.cc | 12 +++++------- grid.h | 4 ++-- vectors.h | 29 +++++++++++++++-------------- 4 files changed, 28 insertions(+), 29 deletions(-) diff --git a/gammapkt.cc b/gammapkt.cc index 1134e30c6..d66b989a5 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -409,8 +409,8 @@ static auto thomson_angle() -> double { return mu; } -[[nodiscard]] static auto scatter_dir(const std::array dir_in, const double cos_theta) - -> std::array +[[nodiscard]] static auto scatter_dir(const std::array dir_in, + const double cos_theta) -> std::array // Routine for scattering a direction through angle theta. { // begin with setting the direction in coordinates where original direction @@ -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, true); + 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, true); + move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist); grid::change_cell(pkt_copy, snext); end_packet = (pkt_copy.type == TYPE_ESCAPE); @@ -1070,7 +1070,7 @@ void guttman_thermalisation(Packet &pkt) { column_densities[i] += grid::get_rho_tmin(mgi) * s_cont; // contribution to the integral } // move packet copy now - move_pkt_withtime(pkt_copy, sdist, true); + move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist); grid::change_cell(pkt_copy, snext); end_packet = (pkt_copy.type == TYPE_ESCAPE); diff --git a/grid.cc b/grid.cc index b0107847c..132be6ed6 100644 --- a/grid.cc +++ b/grid.cc @@ -2462,13 +2462,11 @@ template return -1.; } -static auto get_coordboundary_distances_cylindrical2d(const std::array pkt_pos, - const std::array pkt_dir, - const std::array pktposgridcoord, - const std::array pktvelgridcoord, int cellindex, - const double tstart, const std::array cellcoordmax, - bool neglect_exp) - -> std::tuple, std::array> { +static auto get_coordboundary_distances_cylindrical2d( + const std::array pkt_pos, const std::array pkt_dir, + const std::array pktposgridcoord, const std::array pktvelgridcoord, int cellindex, + const double tstart, const std::array cellcoordmax, + bool neglect_exp) -> std::tuple, std::array> { // 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; diff --git a/grid.h b/grid.h index 95db53000..f09a6f019 100644 --- a/grid.h +++ b/grid.h @@ -130,8 +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 dir, std::array pos, double tstart, int cellindex, - enum cell_boundary *pkt_last_cross, bool neglect_exp = false) - -> std::tuple; + enum cell_boundary *pkt_last_cross, + bool neglect_exp = false) -> std::tuple; [[nodiscard]] inline auto get_elem_abundance(int modelgridindex, int element) -> float // mass fraction of an element (all isotopes combined) diff --git a/vectors.h b/vectors.h index e70be76f1..b21332c49 100644 --- a/vectors.h +++ b/vectors.h @@ -34,15 +34,15 @@ template } template -[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array x, const std::array y) - -> double +[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array x, + const std::array y) -> double // vector dot product { return std::inner_product(x.begin(), x.end(), y.begin(), 0.); } -[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span x, const double t) - -> std::array +[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span x, + const double t) -> std::array // Routine for getting velocity vector of the flow at a position with homologous expansion. { return std::array{x[0] / t, x[1] / t, x[2] / t}; @@ -59,8 +59,8 @@ template return std::array{vec[0] * scalefactor, vec[1] * scalefactor, vec[2] * scalefactor}; } -[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array dir1, const std::array vel) - -> std::array +[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array dir1, + const std::array vel) -> std::array // aberation of angles in special relativity // dir1: direction unit vector in frame1 // vel: velocity of frame2 relative to frame1 @@ -141,17 +141,17 @@ template return doppler_nucmf_on_nurf(dir_rf, get_velocity(pos_rf, prop_time)); } +template constexpr auto move_pkt_withtime(std::span pos_rf, const std::array dir_rf, double &prop_time, const double nu_rf, double &nu_cmf, const double e_rf, double &e_cmf, - const double distance, bool neglect_exp = false) -> double + const double distance) -> double /// Subroutine to move a packet along a straight line (specified by current /// dir vector). The distance moved is in the rest frame. { - double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP; assert_always(distance >= 0); const double nu_cmf_old = nu_cmf; - prop_time += distance / CLIGHT_PROP_DIST_CALC; + prop_time += distance / CLIGHT_ASSUMED; pos_rf[0] += (dir_rf[0] * distance); pos_rf[1] += (dir_rf[1] * distance); @@ -171,9 +171,10 @@ constexpr auto move_pkt_withtime(std::span pos_rf, const std::array 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, - neglect_exp); +template +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); } [[nodiscard]] [[gnu::const]] constexpr auto get_arrive_time(const Packet &pkt) -> double @@ -310,8 +311,8 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance, bool neglec } // Routine to transform the Stokes Parameters from RF to CMF -constexpr auto frame_transform(const std::array n_rf, double *Q, double *U, const std::array v) - -> std::array { +constexpr auto frame_transform(const std::array n_rf, double *Q, double *U, + const std::array v) -> std::array { // Meridian frame in the RF const auto [ref1_rf, ref2_rf] = meridian(n_rf);