Skip to content

Dev/gr pusher axis #126

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

Merged
alisagk merged 3 commits into from
Jul 15, 2025
Merged

Dev/gr pusher axis #126

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
55df3f1
added check for particles close to the axis. the "closeness" is defin…
alisagk Jul 14, 2025
8b692a5
added axis check in the GR deposit
alisagk Jul 14, 2025
a4965f1
minor (deleted useless line)
alisagk Jul 14, 2025
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
25 changes: 13 additions & 12 deletions src/global/utils/numeric.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,18 +80,19 @@ inline constexpr double INV_64 = 0.015625;
#define CROSS_x3(ax1, ax2, ax3, bx1, bx2, bx3) ((ax1) * (bx2) - (ax2) * (bx1))

namespace constant {
inline constexpr std::uint64_t RandomSeed = 0x123456789abcdef0;
inline constexpr double HALF_PI = 1.57079632679489661923;
inline constexpr double PI = 3.14159265358979323846;
inline constexpr double INV_PI = 0.31830988618379067154;
inline constexpr double PI_SQR = 9.86960440108935861882;
inline constexpr double INV_PI_SQR = 0.10132118364233777144;
inline constexpr double TWO_PI = 6.28318530717958647692;
inline constexpr double E = 2.71828182845904523536;
inline constexpr double SQRT2 = 1.41421356237309504880;
inline constexpr double INV_SQRT2 = 0.70710678118654752440;
inline constexpr double SQRT3 = 1.73205080756887729352;
inline constexpr double SMALL_ANGLE = 1e-3;
inline constexpr std::uint64_t RandomSeed = 0x123456789abcdef0;
inline constexpr double HALF_PI = 1.57079632679489661923;
inline constexpr double PI = 3.14159265358979323846;
inline constexpr double INV_PI = 0.31830988618379067154;
inline constexpr double PI_SQR = 9.86960440108935861882;
inline constexpr double INV_PI_SQR = 0.10132118364233777144;
inline constexpr double TWO_PI = 6.28318530717958647692;
inline constexpr double E = 2.71828182845904523536;
inline constexpr double SQRT2 = 1.41421356237309504880;
inline constexpr double INV_SQRT2 = 0.70710678118654752440;
inline constexpr double SQRT3 = 1.73205080756887729352;
inline constexpr double SMALL_ANGLE = 1e-3;
inline constexpr double SMALL_ANGLE_GR = 1e-5;
} // namespace constant

namespace convert {
Expand Down
16 changes: 14 additions & 2 deletions src/kernels/currents_deposit.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,10 +127,22 @@ namespace kernel {
vp);
inv_energy = ONE / math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p)));
} else {
metric.template transform<Idx::D, Idx::U>(xp,
coord_t<Dim::_2D> xp_ { ZERO };
xp_[0] = xp[0];
real_t theta_Cd { xp[1] };
const real_t theta_Ph { metric.template convert<2, Crd::Cd, Crd::Ph>(theta_Cd) };
const real_t small_angle { constant::SMALL_ANGLE_GR };
const auto large_angle { constant::PI - small_angle };
if (theta_Ph < small_angle) {
theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(small_angle);
} else if (theta_Ph >= large_angle) {
theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(large_angle);
}
xp_[1] = theta_Cd;
metric.template transform<Idx::D, Idx::U>(xp_,
{ ux1(p), ux2(p), ux3(p) },
vp);
inv_energy = metric.alpha(xp) / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] +
inv_energy = metric.alpha(xp_) / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] +
ux3(p) * vp[2]);
}
if (Kokkos::isnan(vp[2]) || Kokkos::isinf(vp[2])) {
Expand Down
15 changes: 14 additions & 1 deletion src/kernels/particle_pusher_gr.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -678,6 +678,19 @@ namespace kernel::gr {
xp[0] = i_di_to_Xi(i1(p), dx1(p));
xp[1] = i_di_to_Xi(i2(p), dx2(p));

coord_t<Dim::_2D> xp_ { ZERO };
xp_[0] = xp[0];
real_t theta_Cd { xp[1] };
const real_t theta_Ph { metric.template convert<2, Crd::Cd, Crd::Ph>(theta_Cd) };
const real_t small_angle { constant::SMALL_ANGLE_GR };
const auto large_angle { constant::PI - small_angle };
if (theta_Ph < small_angle) {
theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(small_angle);
} else if (theta_Ph >= large_angle) {
theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(large_angle);
}
xp_[1] = theta_Cd;

vec_t<Dim::_3D> Dp_cntrv { ZERO }, Bp_cntrv { ZERO }, Dp_hat { ZERO },
Bp_hat { ZERO };
interpolateFields(p, Dp_cntrv, Bp_cntrv);
Expand All @@ -694,7 +707,7 @@ namespace kernel::gr {
vp[0] = vp_upd[0];
vp[1] = vp_upd[1];
vp[2] = vp_upd[2];
GeodesicMomentumPush<Massive_t>(Massive_t {}, xp, vp, vp_upd);
GeodesicMomentumPush<Massive_t>(Massive_t {}, xp_, vp, vp_upd);
/* u**_i(n) -> u_i(n + 1/2) */
vp[0] = vp_upd[0];
vp[1] = vp_upd[1];
Expand Down