From a8b6210724e1eb04b12ad0a1b20f8e795df9685f Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 24 Feb 2026 12:18:38 -0500 Subject: [PATCH 01/22] started 4-metric in ks --- src/metrics/kerr_schild.h | 71 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 4df88cf4c..ad915a2fe 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -218,6 +218,77 @@ namespace metric { } } + /** + * metric component with lower indices: g_mu nu + * @param x coordinate array in code units + */ + template + Inline auto g_(const coord_t& x) const -> real_t { + static_assert(i >= 0 && i <= 3, "Invalid index i"); + static_assert(j >= 0 && j <= 3, "Invalid index j"); + if constexpr (i == 1 && j == 1) { + // h_11 + return SQR(dr) * (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)); + } else if constexpr (i == 2 && j == 2) { + // h_22 + return SQR(dtheta) * Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + } else if constexpr (i == 3 && j == 3) { + // h_33 + if constexpr (D == Dim::_2D) { + return A(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * + SQR(math::sin(x[1] * dtheta + x2_min)) / + Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + } else { + return SQR(dphi) * A(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * + SQR(math::sin(x[1] * dtheta + x2_min)) / + Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + } + } else if constexpr ((i == 1 && j == 3) || (i == 3 && j == 1)) { + // h_13 or h_31 + if constexpr (D == Dim::_2D) { + return -dr * a * (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)) * + SQR(math::sin(x[1] * dtheta + x2_min)); + } else { + return -dr * dphi * a * + (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)) * + SQR(math::sin(x[1] * dtheta + x2_min)); + } + } else { + return ZERO; + } + } + + /** + * metric component with upper indices: g^mu nu + * @param x coordinate array in code units + */ + template + Inline auto g(const coord_t& x) const -> real_t { + static_assert(i >= 0 && i <= 3, "Invalid index i"); + static_assert(j >= 0 && j <= 3, "Invalid index j"); + if constexpr (i == 0 && j == 0) { + // g^00 + return - ONE / SQR(alpha(x)); + } else if constexpr ((i == 0 && j == 1) || (i == 1 && j == 0)) { + // g^01 or g^10 + return beta1(x) / SQR(alpha(x)); + } else if constexpr (i == 1 && j == 1) { + // g^11 + return h<1, 1>(x) - SQR(beta1(x)) / SQR(alpha(x)); + } else if constexpr (i == 2 && j == 2) { + // g^22 + return h<2, 2>(x); + } else if constexpr (i == 3 && j == 3) { + // g^33 + return h<3, 3>(x); + } else if constexpr ((i == 1 && j == 3) || (i == 3 && j == 1)) { + // g^13 or g^31 + return h<1, 3>(x); + } else { + return ZERO; + } + } + /** * lapse function * @param x coordinate array in code units From e83ce29359b5d3a383139103a0ada2a46aeb57b0 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 25 Feb 2026 16:37:23 -0500 Subject: [PATCH 02/22] added new functions to support: - calculation of u_0 from u_i componets - 4-metric g_mu nu and g^mu nu - transformation of 4-vectors --- src/metrics/kerr_schild.h | 88 +++++++++++++++++++-------- src/metrics/kerr_schild_0.h | 88 +++++++++++++++++++++++++++ src/metrics/qkerr_schild.h | 115 ++++++++++++++++++++++++++++++++++++ 3 files changed, 265 insertions(+), 26 deletions(-) diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index ad915a2fe..6ea888b77 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -226,33 +226,31 @@ namespace metric { Inline auto g_(const coord_t& x) const -> real_t { static_assert(i >= 0 && i <= 3, "Invalid index i"); static_assert(j >= 0 && j <= 3, "Invalid index j"); - if constexpr (i == 1 && j == 1) { - // h_11 - return SQR(dr) * (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)); - } else if constexpr (i == 2 && j == 2) { - // h_22 - return SQR(dtheta) * Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); - } else if constexpr (i == 3 && j == 3) { - // h_33 + if constexpr (i == 0 && j == 0) { + // g_00 + return - (ONE - z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)); + } else if constexpr ((i == 0 && j == 1) || (i == 1 && j == 0)) { + // g_01 or g_10 + return dr * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + } else if constexpr ((i == 0 && j == 3) || (i == 3 && j == 0)) { + // g_03 or g_30 if constexpr (D == Dim::_2D) { - return A(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * - SQR(math::sin(x[1] * dtheta + x2_min)) / - Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + return - a * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * SQR(math::sin(x[1] * dtheta + x2_min)); } else { - return SQR(dphi) * A(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * - SQR(math::sin(x[1] * dtheta + x2_min)) / - Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + return - dphi_inv * a * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * SQR(math::sin(x[1] * dtheta + x2_min)); } + } else if constexpr (i == 1 && j == 1) { + // g_11 + return SQR(dr) * (1 + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)); + } else if constexpr (i == 2 && j == 2) { + // g_22 + return h_<2, 2>(x); + } else if constexpr (i == 3 && j == 3) { + // g_33 + return h_<3, 3>(x); } else if constexpr ((i == 1 && j == 3) || (i == 3 && j == 1)) { - // h_13 or h_31 - if constexpr (D == Dim::_2D) { - return -dr * a * (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)) * - SQR(math::sin(x[1] * dtheta + x2_min)); - } else { - return -dr * dphi * a * - (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)) * - SQR(math::sin(x[1] * dtheta + x2_min)); - } + // g_13 or g_31 + return h_<1, 3>(x); } else { return ZERO; } @@ -268,13 +266,13 @@ namespace metric { static_assert(j >= 0 && j <= 3, "Invalid index j"); if constexpr (i == 0 && j == 0) { // g^00 - return - ONE / SQR(alpha(x)); + return - (ONE + z(x[0] * dr + x1_min, x[1] * dtheta + x2_min)); } else if constexpr ((i == 0 && j == 1) || (i == 1 && j == 0)) { // g^01 or g^10 - return beta1(x) / SQR(alpha(x)); + return dr_inv * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min); } else if constexpr (i == 1 && j == 1) { // g^11 - return h<1, 1>(x) - SQR(beta1(x)) / SQR(alpha(x)); + return SQR(dr_inv) * Delta(x[0] * dr + x1_min)/ Sigma(x[0] * dr + x1_min, x[1] * dtheta + x2_min); } else if constexpr (i == 2 && j == 2) { // g^22 return h<2, 2>(x); @@ -666,6 +664,44 @@ namespace metric { raise::KernelError(HERE, "Invalid transformation"); } } + + /** + * u_0 from covariant velocity components + */ + Inline auto u_0(const coord_t& xi, const vec_t& u_i) const { + const real_t A { g<0, 0>(xi) }; + const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; + const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; + return (B + math::sqrt(SQR(B) + FOUR * A * C)) / (TWO * A); + } + + /** + * full 4D-vector transformations + */ + template + Inline void transform_4d(const coord_t& xi, + const vec_t& v_in, + vec_t& v_out) const { + static_assert(in != out, "Invalid vector transformation"); + static_assert(in != Idx::XYZ && out != Idx::XYZ, + "Invalid vector transformation: XYZ not allowed in GR"); + if constexpr (in == Idx::D && out == Idx::U) { + // cov -> cntrv + v_out[0] = v_in[0] * g<0, 0>(xi) + v_in[1] * g<0, 1>(xi); + v_out[1] = v_in[0] * g<1, 0>(xi) + v_in[1] * g<1, 1>(xi) + v_in[3] * g<1, 3>(xi); + v_out[2] = v_in[2] * g<2, 2>(xi); + v_out[3] = v_in[1] * g<3, 1>(xi) + v_in[3] * g<3, 3>(xi); + } else if constexpr (in == Idx::U && out == Idx::D) { + // cntrv -> cov + v_out[0] = v_in[0] * g_<0, 0>(xi) + v_in[1] * g_<0, 1>(xi) + v_in[3] * g_<0, 3>(xi); + v_out[1] = v_in[0] * g_<1, 0>(xi) + v_in[1] * g_<1, 1>(xi) + v_in[3] * g_<1, 3>(xi); + v_out[2] = v_in[2] * g_<2, 2>(xi); + v_out[3] = v_in[0] * g_<3, 0>(xi) + v_in[1] * g_<3, 1>(xi) + v_in[3] * g_<3, 3>(xi); + } else { + raise::KernelError(HERE, "Invalid transformation"); + } + } }; } // namespace metric diff --git a/src/metrics/kerr_schild_0.h b/src/metrics/kerr_schild_0.h index 66834a4c7..7e0af95ee 100644 --- a/src/metrics/kerr_schild_0.h +++ b/src/metrics/kerr_schild_0.h @@ -177,6 +177,56 @@ namespace metric { } } + /** + * metric component with lower indices: g_mu nu + * @param x coordinate array in code units + */ + template + Inline auto g_(const coord_t& x) const -> real_t { + static_assert(i >= 0 && i <= 3, "Invalid index i"); + static_assert(j >= 0 && j <= 3, "Invalid index j"); + if constexpr (i == 0 && j == 0) { + // g_00 + return - ONE; + } else if constexpr (i == 1 && j == 1) { + // g_11 + return h_<1, 1>(x);; + } else if constexpr (i == 2 && j == 2) { + // g_22 + return h_<2, 2>(x); + } else if constexpr (i == 3 && j == 3) { + // g_33 + return h_<3, 3>(x); + } else { + return ZERO; + } + } + + /** + * metric component with upper indices: g^mu nu + * @param x coordinate array in code units + */ + template + Inline auto g(const coord_t& x) const -> real_t { + static_assert(i >= 0 && i <= 3, "Invalid index i"); + static_assert(j >= 0 && j <= 3, "Invalid index j"); + if constexpr (i == 0 && j == 0) { + // g^00 + return - ONE; + } else if constexpr (i == 1 && j == 1) { + // g^11 + return h<1, 1>(x); + } else if constexpr (i == 2 && j == 2) { + // g^22 + return h<2, 2>(x); + } else if constexpr (i == 3 && j == 3) { + // g^33 + return h<3, 3>(x); + } else { + return ZERO; + } + } + /** * lapse function * @param x coordinate array in code units @@ -457,6 +507,44 @@ namespace metric { raise::KernelError(HERE, "Invalid transformation"); } } + + /** + * u_0 from covariant velocity components + */ + Inline auto u_0(const coord_t& xi, const vec_t& u_i) const { + const real_t A { g<0, 0>(xi) }; + const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; + const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; + return (B + math::sqrt(SQR(B) + FOUR * A * C)) / (TWO * A); + } + + /** + * full 4D-vector transformations + */ + template + Inline void transform_4d(const coord_t& xi, + const vec_t& v_in, + vec_t& v_out) const { + static_assert(in != out, "Invalid vector transformation"); + static_assert(in != Idx::XYZ && out != Idx::XYZ, + "Invalid vector transformation: XYZ not allowed in GR"); + if constexpr (in == Idx::D && out == Idx::U) { + // cov -> cntrv + v_out[0] = v_in[0] * g<0, 0>(xi) + v_in[1] * g<0, 1>(xi); + v_out[1] = v_in[0] * g<1, 0>(xi) + v_in[1] * g<1, 1>(xi) + v_in[3] * g<1, 3>(xi); + v_out[2] = v_in[2] * g<2, 2>(xi); + v_out[3] = v_in[1] * g<3, 1>(xi) + v_in[3] * g<3, 3>(xi); + } else if constexpr (in == Idx::U && out == Idx::D) { + // cntrv -> cov + v_out[0] = v_in[0] * g_<0, 0>(xi) + v_in[1] * g_<0, 1>(xi) + v_in[3] * g_<0, 3>(xi); + v_out[1] = v_in[0] * g_<1, 0>(xi) + v_in[1] * g_<1, 1>(xi) + v_in[3] * g_<1, 3>(xi); + v_out[2] = v_in[2] * g_<2, 2>(xi); + v_out[3] = v_in[0] * g_<3, 0>(xi) + v_in[1] * g_<3, 1>(xi) + v_in[3] * g_<3, 3>(xi); + } else { + raise::KernelError(HERE, "Invalid transformation"); + } + } }; } // namespace metric diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index 90f366efb..4af96d4c7 100644 --- a/src/metrics/qkerr_schild.h +++ b/src/metrics/qkerr_schild.h @@ -235,6 +235,82 @@ namespace metric { } } + /** + * metric component with lower indices: g_mu nu + * @param x coordinate array in code units + */ + template + Inline auto g_(const coord_t& x) const -> real_t { + static_assert(i >= 0 && i <= 3, "Invalid index i"); + static_assert(j >= 0 && j <= 3, "Invalid index j"); + if constexpr (i == 0 && j == 0) { + // g_00 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + return - (ONE - z(r0 + math::exp(x[0] * dchi + chi_min), theta)); + } else if constexpr ((i == 0 && j == 1) || (i == 1 && j == 0)) { + // g_01 or g_10 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + return dchi * math::exp(x[0] * dchi + chi_min) * z(r0 + math::exp(x[0] * dchi + chi_min), theta); + } else if constexpr ((i == 0 && j == 3) || (i == 3 && j == 0)) { + // g_03 or g_30 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + if constexpr (D == Dim::_2D) { + return - a * z(r0 + math::exp(x[0] * dchi + chi_min)) * SQR(math::sin(theta)); + } else { + return - dphi_inv * a * z(r0 + math::exp(x[0] * dchi + chi_min)) * SQR(math::sin(theta)); + } + } else if constexpr (i == 1 && j == 1) { + // g_11 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + return SQR(dchi) * math::exp(TWO * (x[0] * dchi + chi_min)) * (1 + z(r0 + math::exp(x[0] * dchi + chi_min), theta)); + } else if constexpr (i == 2 && j == 2) { + // g_22 + return h_<2, 2>(x); + } else if constexpr (i == 3 && j == 3) { + // g_33 + return h_<3, 3>(x); + } else if constexpr ((i == 1 && j == 3) || (i == 3 && j == 1)) { + // g_13 or g_31 + return h_<1, 3>(x); + } else { + return ZERO; + } + } + + /** + * metric component with upper indices: g^mu nu + * @param x coordinate array in code units + */ + template + Inline auto g(const coord_t& x) const -> real_t { + static_assert(i >= 0 && i <= 3, "Invalid index i"); + static_assert(j >= 0 && j <= 3, "Invalid index j"); + if constexpr (i == 0 && j == 0) { + // g^00 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + return - (ONE + z(r0 + math::exp(x[0] * dchi + chi_min), theta)); + } else if constexpr ((i == 0 && j == 1) || (i == 1 && j == 0)) { + // g^01 or g^10 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + return (math::exp(-(x[0] * dchi + chi_min)) * dchi_inv) * z(r0 + math::exp(x[0] * dchi + chi_min), theta); + } else if constexpr (i == 1 && j == 1) { + // g^11 + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + return (math::exp(-TWO * (x[0] * dchi + chi_min)) / SQR(dchi)) * Delta(r0 + math::exp(x[0] * dchi + chi_min))/ Sigma(r0 + math::exp(x[0] * dchi + chi_min), theta); + } else if constexpr (i == 2 && j == 2) { + // g^22 + return h<2, 2>(x); + } else if constexpr (i == 3 && j == 3) { + // g^33 + return h<3, 3>(x); + } else if constexpr ((i == 1 && j == 3) || (i == 3 && j == 1)) { + // g^13 or g^31 + return h<1, 3>(x); + } else { + return ZERO; + } + } + /** * lapse function * @param x coordinate array in code units @@ -663,6 +739,45 @@ namespace metric { } } + + /** + * u_0 from covariant velocity components + */ + Inline auto u_0(const coord_t& xi, const vec_t& u_i) const { + const real_t A { g<0, 0>(xi) }; + const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; + const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; + return (B + math::sqrt(SQR(B) + FOUR * A * C)) / (TWO * A); + } + + /** + * full 4D-vector transformations + */ + template + Inline void transform_4d(const coord_t& xi, + const vec_t& v_in, + vec_t& v_out) const { + static_assert(in != out, "Invalid vector transformation"); + static_assert(in != Idx::XYZ && out != Idx::XYZ, + "Invalid vector transformation: XYZ not allowed in GR"); + if constexpr (in == Idx::D && out == Idx::U) { + // cov -> cntrv + v_out[0] = v_in[0] * g<0, 0>(xi) + v_in[1] * g<0, 1>(xi); + v_out[1] = v_in[0] * g<1, 0>(xi) + v_in[1] * g<1, 1>(xi) + v_in[3] * g<1, 3>(xi); + v_out[2] = v_in[2] * g<2, 2>(xi); + v_out[3] = v_in[1] * g<3, 1>(xi) + v_in[3] * g<3, 3>(xi); + } else if constexpr (in == Idx::U && out == Idx::D) { + // cntrv -> cov + v_out[0] = v_in[0] * g_<0, 0>(xi) + v_in[1] * g_<0, 1>(xi) + v_in[3] * g_<0, 3>(xi); + v_out[1] = v_in[0] * g_<1, 0>(xi) + v_in[1] * g_<1, 1>(xi) + v_in[3] * g_<1, 3>(xi); + v_out[2] = v_in[2] * g_<2, 2>(xi); + v_out[3] = v_in[0] * g_<3, 0>(xi) + v_in[1] * g_<3, 1>(xi) + v_in[3] * g_<3, 3>(xi); + } else { + raise::KernelError(HERE, "Invalid transformation"); + } + } + private: /* Specific for Quasi-Spherical metric with angle stretching ------------ */ /** From 204d37ab4938bdb962077d72abc5e9a0eaf9f0c7 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 25 Feb 2026 17:03:16 -0500 Subject: [PATCH 03/22] bugs in 4-metric implemetations --- src/metrics/kerr_schild.h | 4 ++-- src/metrics/kerr_schild_0.h | 2 +- src/metrics/qkerr_schild.h | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 6ea888b77..5e0dc9694 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -237,7 +237,7 @@ namespace metric { if constexpr (D == Dim::_2D) { return - a * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * SQR(math::sin(x[1] * dtheta + x2_min)); } else { - return - dphi_inv * a * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * SQR(math::sin(x[1] * dtheta + x2_min)); + return - dphi * a * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * SQR(math::sin(x[1] * dtheta + x2_min)); } } else if constexpr (i == 1 && j == 1) { // g_11 @@ -673,7 +673,7 @@ namespace metric { const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; - return (B + math::sqrt(SQR(B) + FOUR * A * C)) / (TWO * A); + return (B + math::sqrt(SQR(B) - FOUR * A * C)) / (TWO * A); } /** diff --git a/src/metrics/kerr_schild_0.h b/src/metrics/kerr_schild_0.h index 7e0af95ee..2b8510fa1 100644 --- a/src/metrics/kerr_schild_0.h +++ b/src/metrics/kerr_schild_0.h @@ -516,7 +516,7 @@ namespace metric { const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; - return (B + math::sqrt(SQR(B) + FOUR * A * C)) / (TWO * A); + return (B + math::sqrt(SQR(B) - FOUR * A * C)) / (TWO * A); } /** diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index 4af96d4c7..0f4615bf2 100644 --- a/src/metrics/qkerr_schild.h +++ b/src/metrics/qkerr_schild.h @@ -255,9 +255,9 @@ namespace metric { // g_03 or g_30 const real_t theta { eta2theta(x[1] * deta + eta_min) }; if constexpr (D == Dim::_2D) { - return - a * z(r0 + math::exp(x[0] * dchi + chi_min)) * SQR(math::sin(theta)); + return - a * z(r0 + math::exp(x[0] * dchi + chi_min), theta) * SQR(math::sin(theta)); } else { - return - dphi_inv * a * z(r0 + math::exp(x[0] * dchi + chi_min)) * SQR(math::sin(theta)); + return - dphi * a * z(r0 + math::exp(x[0] * dchi + chi_min)) * SQR(math::sin(theta)); } } else if constexpr (i == 1 && j == 1) { // g_11 @@ -748,7 +748,7 @@ namespace metric { const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; - return (B + math::sqrt(SQR(B) + FOUR * A * C)) / (TWO * A); + return (B + math::sqrt(SQR(B) - FOUR * A * C)) / (TWO * A); } /** From d9de68e5954e7efcf4319803dcfceb88bd4e5c59 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 25 Feb 2026 17:04:35 -0500 Subject: [PATCH 04/22] bug in qks g_ implementation --- src/metrics/qkerr_schild.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index 0f4615bf2..1636ca737 100644 --- a/src/metrics/qkerr_schild.h +++ b/src/metrics/qkerr_schild.h @@ -257,7 +257,7 @@ namespace metric { if constexpr (D == Dim::_2D) { return - a * z(r0 + math::exp(x[0] * dchi + chi_min), theta) * SQR(math::sin(theta)); } else { - return - dphi * a * z(r0 + math::exp(x[0] * dchi + chi_min)) * SQR(math::sin(theta)); + return - dphi * a * z(r0 + math::exp(x[0] * dchi + chi_min), theta) * SQR(math::sin(theta)); } } else if constexpr (i == 1 && j == 1) { // g_11 From ae8bd688ef2c240689dc68b66368a88b59bdf3db Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 12:58:31 -0500 Subject: [PATCH 05/22] option for u_0 calculation for photons and massive particles --- src/metrics/kerr_schild.h | 4 ++-- src/metrics/kerr_schild_0.h | 4 ++-- src/metrics/qkerr_schild.h | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 5e0dc9694..7670549b1 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -668,11 +668,11 @@ namespace metric { /** * u_0 from covariant velocity components */ - Inline auto u_0(const coord_t& xi, const vec_t& u_i) const { + Inline auto u_0(const coord_t& xi, const vec_t& u_i, const real_t norm) const { const real_t A { g<0, 0>(xi) }; const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + - g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + norm }; return (B + math::sqrt(SQR(B) - FOUR * A * C)) / (TWO * A); } diff --git a/src/metrics/kerr_schild_0.h b/src/metrics/kerr_schild_0.h index 2b8510fa1..e1debf689 100644 --- a/src/metrics/kerr_schild_0.h +++ b/src/metrics/kerr_schild_0.h @@ -511,11 +511,11 @@ namespace metric { /** * u_0 from covariant velocity components */ - Inline auto u_0(const coord_t& xi, const vec_t& u_i) const { + Inline auto u_0(const coord_t& xi, const vec_t& u_i, const real_t norm) const { const real_t A { g<0, 0>(xi) }; const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + - g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + norm }; return (B + math::sqrt(SQR(B) - FOUR * A * C)) / (TWO * A); } diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index 1636ca737..beab85b0a 100644 --- a/src/metrics/qkerr_schild.h +++ b/src/metrics/qkerr_schild.h @@ -743,11 +743,11 @@ namespace metric { /** * u_0 from covariant velocity components */ - Inline auto u_0(const coord_t& xi, const vec_t& u_i) const { + Inline auto u_0(const coord_t& xi, const vec_t& u_i, const real_t norm) const { const real_t A { g<0, 0>(xi) }; const real_t B { - TWO * g<0, 1>(xi) * u_i[0] }; const real_t C { g<1, 1>(xi) * SQR(u_i[0]) + g<2, 2>(xi) * SQR(u_i[1]) + - g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + ONE }; + g<3, 3>(xi) * SQR(u_i[2]) + TWO * g<1, 3>(xi) * u_i[0] * u_i[2] + norm }; return (B + math::sqrt(SQR(B) - FOUR * A * C)) / (TWO * A); } From 1f81a5503030b59e01b025240b184c06be434fd2 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 13:35:55 -0500 Subject: [PATCH 06/22] removed dublicate computation of FldsID::V --- src/kernels/particle_moments.hpp | 59 -------------------------------- 1 file changed, 59 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 023eb55a5..be825409d 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -235,65 +235,6 @@ namespace kernel { // for other cases, use the `contrib` defined above coeff = contrib; } - - if constexpr (F == FldsID::V) { - real_t gamma { ZERO }; - // for stress-energy tensor - vec_t u_Phys { ZERO }; - if constexpr (S == SimEngine::SRPIC) { - // SR - // stress-energy tensor for SR is computed in the tetrad (hatted) basis - if constexpr (M::CoordType == Coord::Cart) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); - } else { - static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D"); - coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); - if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); - } else { - x_Code[2] = phi(p); - } - metric.template transform_xyz( - x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); - } - if (mass == ZERO) { - gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); - } else { - gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); - } - } else { - // GR - // stress-energy tensor for GR is computed in contravariant basis - static_assert(D != Dim::_1D, "GRPIC 1D"); - coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); - if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); - } - vec_t u_Cntrv { ZERO }; - // compute u_i u^i for energy - metric.template transform(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Cntrv); - gamma = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + u_Cntrv[2] * ux3(p); - if (mass == ZERO) { - gamma = math::sqrt(gamma); - } else { - gamma = math::sqrt(ONE + gamma); - } - metric.template transform(x_Code, u_Cntrv, u_Phys); - } - // compute the corresponding moment - coeff = u_Phys[c1 - 1] / gamma; - } - if constexpr (F != FldsID::Nppc) { // for nppc calculation ... // ... do not take volume, weights or smoothing into account From 43fbe341122f71221cd7efbea5852fbf77560794 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 13:56:24 -0500 Subject: [PATCH 07/22] fixes bugs in tensor calculation: 1. all tensor components should have mass (T^00 did not have it) 2. all tensor components should have only one mass (T^ij had m^2) Fixed by multiplying the coefficient by mass at the beginning only --- src/kernels/particle_moments.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index be825409d..f9554e38a 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -194,13 +194,14 @@ namespace kernel { metric.template transform(x_Code, u_Cntrv, u_Phys); } // compute the corresponding moment - coeff = ONE / energy; + // T^μν = mass * u^μ * u^ν / u^0 + coeff = ((mass == ZERO) ? ONE : mass) / energy; #pragma unroll for (const auto& c : { c1, c2 }) { if (c == 0) { - coeff *= energy; + coeff *= energy; // multiply by u^0 } else { - coeff *= ((mass == ZERO) ? ONE : mass) * u_Phys[c - 1]; + coeff *= u_Phys[c - 1]; // multiply by u^c } } } else if constexpr (F == FldsID::V) { From 88e0f48a7cf188244609eaa2eb604217a53d31fb Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 14:03:59 -0500 Subject: [PATCH 08/22] removed unnecessary commets --- src/kernels/particle_moments.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index f9554e38a..984634c98 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -199,9 +199,9 @@ namespace kernel { #pragma unroll for (const auto& c : { c1, c2 }) { if (c == 0) { - coeff *= energy; // multiply by u^0 + coeff *= energy; } else { - coeff *= u_Phys[c - 1]; // multiply by u^c + coeff *= u_Phys[c - 1]; } } } else if constexpr (F == FldsID::V) { From e76c03dff0caa302bf46d1fdc102a8aa59050a09 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 14:39:27 -0500 Subject: [PATCH 09/22] stress energy tensor now uses the correct contravariant velocity --- src/kernels/particle_moments.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 984634c98..7e53dc720 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -191,6 +191,18 @@ namespace kernel { } else { energy = mass * math::sqrt(ONE + energy); } + // raise full covariant 4-vector to get correct contravariant u^i + // u^i != h^{ij} u_j + const real_t norm { (mass == ZERO) ? ZERO : ONE }; + const real_t u_0_cov { metric.u_0(x_Code, { ux1(p), ux2(p), ux3(p) }, norm) }; + vec_t u_cntrv_4d { ZERO }; + metric.template transform_4d( + x_Code, + { u_0_cov, ux1(p), ux2(p), ux3(p) }, + u_cntrv_4d); + u_Cntrv[0] = u_cntrv_4d[1]; + u_Cntrv[1] = u_cntrv_4d[2]; + u_Cntrv[2] = u_cntrv_4d[3]; metric.template transform(x_Code, u_Cntrv, u_Phys); } // compute the corresponding moment From 20a3a05eaffb755b25f8647af346f4a5b999d8d8 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 15:58:29 -0500 Subject: [PATCH 10/22] fixed m^2 and simplified calculation --- src/kernels/particle_moments.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 7e53dc720..cc772a000 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -206,14 +206,14 @@ namespace kernel { metric.template transform(x_Code, u_Cntrv, u_Phys); } // compute the corresponding moment - // T^μν = mass * u^μ * u^ν / u^0 - coeff = ((mass == ZERO) ? ONE : mass) / energy; + // T^μν = energy * v^{c1} * v^{c2}, + // where v^0 = 1, v^i = u^i / u^0 (coordinate 3-velocity) + const real_t inv_u0 { (mass == ZERO) ? ONE : mass / energy }; + coeff = energy; #pragma unroll for (const auto& c : { c1, c2 }) { - if (c == 0) { - coeff *= energy; - } else { - coeff *= u_Phys[c - 1]; + if (c != 0) { + coeff *= u_Phys[c - 1] * inv_u0; } } } else if constexpr (F == FldsID::V) { From 93c3722d9ce74b9f0119825c3e7b9920f9cd6404 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 16:20:14 -0500 Subject: [PATCH 11/22] added correct u^0 in GR; cut energy definition out -- now uses only u^0 --- src/kernels/particle_moments.hpp | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index cc772a000..870bb19e9 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -139,7 +139,7 @@ namespace kernel { } real_t coeff { ZERO }; if constexpr (F == FldsID::T) { - real_t energy { ZERO }; + real_t u0 { ZERO }; // for stress-energy tensor vec_t u_Phys { ZERO }; if constexpr (S == SimEngine::SRPIC) { @@ -165,10 +165,9 @@ namespace kernel { u_Phys); } if (mass == ZERO) { - energy = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); + u0 = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); } else { - energy = mass * - math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); + u0 = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); } } else { // GR @@ -181,16 +180,6 @@ namespace kernel { x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); } vec_t u_Cntrv { ZERO }; - // compute u_i u^i for energy - metric.template transform(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Cntrv); - energy = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + u_Cntrv[2] * ux3(p); - if (mass == ZERO) { - energy = math::sqrt(energy); - } else { - energy = mass * math::sqrt(ONE + energy); - } // raise full covariant 4-vector to get correct contravariant u^i // u^i != h^{ij} u_j const real_t norm { (mass == ZERO) ? ZERO : ONE }; @@ -200,20 +189,21 @@ namespace kernel { x_Code, { u_0_cov, ux1(p), ux2(p), ux3(p) }, u_cntrv_4d); + // in GR (3+1): u^0 = Gamma/alpha (includes lapse) + u0 = u_cntrv_4d[0]; u_Cntrv[0] = u_cntrv_4d[1]; u_Cntrv[1] = u_cntrv_4d[2]; u_Cntrv[2] = u_cntrv_4d[3]; metric.template transform(x_Code, u_Cntrv, u_Phys); } // compute the corresponding moment - // T^μν = energy * v^{c1} * v^{c2}, + // T^μν = m * u^0 * v^{c1} * v^{c2}, // where v^0 = 1, v^i = u^i / u^0 (coordinate 3-velocity) - const real_t inv_u0 { (mass == ZERO) ? ONE : mass / energy }; - coeff = energy; + coeff = ((mass == ZERO) ? ONE : mass) * u0; #pragma unroll for (const auto& c : { c1, c2 }) { if (c != 0) { - coeff *= u_Phys[c - 1] * inv_u0; + coeff *= u_Phys[c - 1] / u0; } } } else if constexpr (F == FldsID::V) { From fdb99f050e63d2d4005ab3445f807b991874206b Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Tue, 3 Mar 2026 16:25:52 -0500 Subject: [PATCH 12/22] cleaned some variables --- src/kernels/particle_moments.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 870bb19e9..80c6f5f12 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -179,22 +179,22 @@ namespace kernel { if constexpr (D == Dim::_3D) { x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); } - vec_t u_Cntrv { ZERO }; // raise full covariant 4-vector to get correct contravariant u^i // u^i != h^{ij} u_j - const real_t norm { (mass == ZERO) ? ZERO : ONE }; - const real_t u_0_cov { metric.u_0(x_Code, { ux1(p), ux2(p), ux3(p) }, norm) }; + const real_t u_0_cov { metric.u_0(x_Code, + { ux1(p), ux2(p), ux3(p) }, + (mass == ZERO) ? ZERO : ONE) }; vec_t u_cntrv_4d { ZERO }; metric.template transform_4d( x_Code, { u_0_cov, ux1(p), ux2(p), ux3(p) }, u_cntrv_4d); - // in GR (3+1): u^0 = Gamma/alpha (includes lapse) - u0 = u_cntrv_4d[0]; - u_Cntrv[0] = u_cntrv_4d[1]; - u_Cntrv[1] = u_cntrv_4d[2]; - u_Cntrv[2] = u_cntrv_4d[3]; - metric.template transform(x_Code, u_Cntrv, u_Phys); + // in GR: u^0 = Gamma/alpha + u0 = u_cntrv_4d[0]; + metric.template transform( + x_Code, + { u_cntrv_4d[1], u_cntrv_4d[2], u_cntrv_4d[3] }, + u_Phys); } // compute the corresponding moment // T^μν = m * u^0 * v^{c1} * v^{c2}, From fa6bfd352dbc94a566d9c69b24a9508efe4ac300 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 4 Mar 2026 13:26:49 -0500 Subject: [PATCH 13/22] Eckart frame transformation --- src/framework/domain/output.cpp | 75 ++++++++- src/global/global.h | 3 +- src/kernels/particle_moments.hpp | 252 ++++++++++++++++++++++++++++--- src/output/fields.cpp | 33 +++- 4 files changed, 328 insertions(+), 35 deletions(-) diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 1d29c4102..9601b1ddb 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -487,7 +487,16 @@ namespace ntt { local_domain->fields.bckp, c); } else if (fld.id() == FldsID::V) { - if constexpr (S != SimEngine::GRPIC) { + if constexpr (S == SimEngine::SRPIC) { + ComputeMoments(params, + local_domain->mesh, + local_domain->species, + fld.species, + fld.comp[0], + local_domain->fields.bckp, + c); + } else if constexpr (S == SimEngine::GRPIC) { + // Single component of 4-velocity (user requested e.g., "v_0") ComputeMoments(params, local_domain->mesh, local_domain->species, @@ -495,8 +504,6 @@ namespace ntt { fld.comp[0], local_domain->fields.bckp, c); - } else { - raise::Error("Bulk velocity not supported for GRPIC", HERE); } } else { raise::Error("Wrong moment requested for output", HERE); @@ -577,8 +584,15 @@ namespace ntt { fld.comp[i], local_domain->fields.bckp, c); - } else { - raise::Error("Bulk velocity not supported for GRPIC", HERE); + } else if constexpr (S == SimEngine::GRPIC) { + // Single component from 3-vector spec (legacy), deposit as part of 4-vector + ComputeMoments(params, + local_domain->mesh, + local_domain->species, + fld.species, + fld.comp[i], + local_domain->fields.bckp, + c); } } else { raise::Error("Wrong moment requested for output", HERE); @@ -678,6 +692,57 @@ namespace ntt { local_domain->mesh.metric)); } } + } else if (fld.comp.size() == 4) { // 4-vector + if constexpr (S == SimEngine::GRPIC) { + if (fld.is_moment() && fld.id() == FldsID::V) { + // Compute 4-velocity: V^μ (u^0, u^1, u^2, u^3) + for (auto i = 0; i < 4; ++i) { + names.push_back(fld.name(i)); + addresses.push_back(i); + const auto c = static_cast(addresses[i]); + raise::ErrorIf(fld.comp[i].size() != 1, + "Wrong # of components requested for 4-velocity", + HERE); + ComputeMoments(params, + local_domain->mesh, + local_domain->species, + fld.species, + fld.comp[i], + local_domain->fields.bckp, + c); + } + // Synchronize all 4 components + SynchronizeFields(*local_domain, Comm::Bckp, { addresses[0], addresses[3] + 1 }); + // Normalize 4-momentum flux: V^μ = N^μ / sqrt(-N_ν N^ν) + // (computed in coordinate contravariant basis) + Kokkos::parallel_for("Normalize4VelocityByNorm", + local_domain->mesh.rangeActiveCells(), + kernel::Normalize4VelocityByNorm_kernel( + local_domain->fields.bckp, + local_domain->fields.bckp, + addresses[0], // c_u0 (column for u^0) + addresses[1], // c_u1 (column for u^1) + addresses[2], // c_u2 (column for u^2) + addresses[3], // c_u3 (column for u^3) + local_domain->mesh.metric)); + SynchronizeFields(*local_domain, Comm::Bckp, { addresses[0], addresses[3] + 1 }); + // Transform spatial components to physical basis for output + // u^0 (Gamma/alpha) remains unitless, only u^i transform + Kokkos::parallel_for("Transform4VelocitySpatialToPhysical", + local_domain->mesh.rangeActiveCells(), + kernel::Transform4VelocitySpatialToPhysical_kernel( + local_domain->fields.bckp, + addresses[1], // c_u1 (column for u^1) + addresses[2], // c_u2 (column for u^2) + addresses[3], // c_u3 (column for u^3) + local_domain->mesh.metric)); + } else { + raise::Error("4-vector output only supported for V (bulk velocity) moment in GRPIC", + HERE); + } + } else { + raise::Error("4-vector output only supported for GRPIC", HERE); + } } else if (fld.comp.size() == 6) { // tensor raise::ErrorIf(not fld.is_moment() or fld.id() != FldsID::T, "Only T tensor has 6 components", diff --git a/src/global/global.h b/src/global/global.h index 45da4f043..d6edfb467 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -168,7 +168,8 @@ namespace ntt { enum Dimension : unsigned short { _1D = 1, _2D = 2, - _3D = 3 + _3D = 3, + _4D = 4 }; enum class CellLayer { diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 80c6f5f12..1dc59d6af 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -43,8 +43,6 @@ namespace kernel { class ParticleMoments_kernel { static constexpr auto D = M::Dim; - static_assert(!((S == SimEngine::GRPIC) && (F == FldsID::V)), - "Bulk velocity not supported for GRPIC"); static_assert((F == FldsID::Rho) || (F == FldsID::Charge) || (F == FldsID::N) || (F == FldsID::Nppc) || (F == FldsID::T) || (F == FldsID::V), "Invalid field ID"); @@ -207,33 +205,71 @@ namespace kernel { } } } else if constexpr (F == FldsID::V) { - real_t gamma { ZERO }; - // for bulk 3vel (tetrad basis) - vec_t u_Phys { ZERO }; - if constexpr (M::CoordType == Coord::Cart) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); + if constexpr (S == SimEngine::SRPIC) { + // SR: bulk velocity + real_t gamma { ZERO }; + // for bulk 3vel (tetrad basis) + vec_t u_Phys { ZERO }; + if constexpr (M::CoordType == Coord::Cart) { + u_Phys[0] = ux1(p); + u_Phys[1] = ux2(p); + u_Phys[2] = ux3(p); + } else { + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); + x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + if constexpr (D == Dim::_3D) { + x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + } else { + x_Code[2] = phi(p); + } + metric.template transform_xyz(x_Code, + { ux1(p), ux2(p), ux3(p) }, + u_Phys); + } + if (mass == ZERO) { + gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); + } else { + gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); + } + // compute the corresponding moment + coeff = (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / gamma; } else { - coord_t x_Code { ZERO }; + // GR: Eckart frame flux N^μ = m * u^μ / u^0 + static_assert(D != Dim::_1D, "GRPIC 1D"); + vec_t u_Phys { ZERO }; + real_t u0 { ZERO }; + coord_t x_Code { ZERO }; x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); if constexpr (D == Dim::_3D) { x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + } + // raise full covariant 4-vector to get correct contravariant u^i + // u^i != h^{ij} u_j + const real_t u_0_cov { metric.u_0(x_Code, + { ux1(p), ux2(p), ux3(p) }, + (mass == ZERO) ? ZERO : ONE) }; + vec_t u_cntrv_4d { ZERO }; + metric.template transform_4d( + x_Code, + { u_0_cov, ux1(p), ux2(p), ux3(p) }, + u_cntrv_4d); + // in GR: u^0 = Gamma/alpha + u0 = u_cntrv_4d[0]; + metric.template transform( + x_Code, + { u_cntrv_4d[1], u_cntrv_4d[2], u_cntrv_4d[3] }, + u_Phys); + // Deposit flux N^μ = mass * u^μ / u^0 + if (c1 == 0) { + // u^0 component + coeff = (mass == ZERO ? ONE : mass); } else { - x_Code[2] = phi(p); + // u^i component: mass * u^i / u^0 + coeff = (mass == ZERO ? ONE : mass) * u_cntrv_4d[c1] / u0; } - metric.template transform_xyz(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); - } - if (mass == ZERO) { - gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); - } else { - gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); } - // compute the corresponding moment - coeff = (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / gamma; } else { // for other cases, use the `contrib` defined above coeff = contrib; @@ -396,6 +432,180 @@ namespace kernel { } }; + template + requires metric::traits::HasD && metric::traits::HasTransform + class Normalize4VelocityByNorm_kernel { + // Normalizes 4-momentum flux to Eckart frame velocity + // V^μ = N^μ / sqrt(-N_ν N^ν) + const ndfield_t Flux; // momentum flux N^μ + ndfield_t Vector; // Eckart 4-velocity + const unsigned short c_u0, c_u1, c_u2, c_u3; // 4-velocity component indices + const M metric; + + public: + Normalize4VelocityByNorm_kernel(const ndfield_t& flux, + const ndfield_t& vector, + unsigned short cu0, + unsigned short cu1, + unsigned short cu2, + unsigned short cu3, + const M& metric) + : Flux { flux } + , Vector { vector } + , c_u0 { cu0 } + , c_u1 { cu1 } + , c_u2 { cu2 } + , c_u3 { cu3 } + , metric { metric } { + raise::ErrorIf(c_u0 >= N or c_u1 >= N or c_u2 >= N or c_u3 >= N, + "Invalid component index", + HERE); + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1) + HALF; + x_Code[1] = static_cast(i2) + HALF; + + vec_t N_cntrv { ZERO }; + N_cntrv[0] = Flux(i1, i2, c_u0); + N_cntrv[1] = Flux(i1, i2, c_u1); + N_cntrv[2] = Flux(i1, i2, c_u2); + N_cntrv[3] = Flux(i1, i2, c_u3); + + vec_t N_cov { ZERO }; + // Compute N_i = g_ij N^j + metric.template transform_4d(x_Code, N_cntrv, N_cov); + + // Compute N_ν N^ν = -(N^0)^2 + g_ij N^i N^j + real_t N_norm_sq { N_cov[0] * N_cntrv[0] + N_cov[1] * N_cntrv[1] + N_cov[2] * N_cntrv[2] + N_cov[3] * N_cntrv[3] }; + real_t norm = math::sqrt(-N_norm_sq); + + if (not cmp::AlmostZero(norm)) { + Vector(i1, i2, c_u0) = N_cntrv[0] / norm; + Vector(i1, i2, c_u1) = N_cntrv[1] / norm; + Vector(i1, i2, c_u2) = N_cntrv[2] / norm; + Vector(i1, i2, c_u3) = N_cntrv[3] / norm; + } + } else { + raise::KernelError( + HERE, + "2D implementation of Normalize4VelocityByNorm_kernel called for non-2D"); + } + } + + Inline void operator()(index_t i1, index_t i2, index_t i3) const { + if constexpr (D == Dim::_3D) { + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1) + HALF; + x_Code[1] = static_cast(i2) + HALF; + x_Code[2] = static_cast(i3) + HALF; + + vec_t N_cntrv { ZERO }; + N_cntrv[0] = Flux(i1, i2, i3, c_u0); // N^0 + N_cntrv[1] = Flux(i1, i2, i3, c_u1); // N^1 + N_cntrv[2] = Flux(i1, i2, i3, c_u2); // N^2 + N_cntrv[3] = Flux(i1, i2, i3, c_u3); // N^3 + + vec_t N_cov { ZERO }; + // Compute N_i = g_ij N^j + metric.template transform_4d(x_Code, N_cntrv, N_cov); + + // Compute N_ν N^ν = -(N^0)^2 + g_ij N^i N^j + real_t N_norm_sq { N_cov[0] * N_cntrv[0] + N_cov[1] * N_cntrv[1] + N_cov[2] * N_cntrv[2] + N_cov[3] * N_cntrv[3] }; + + real_t norm = math::sqrt(-N_norm_sq); // sqrt(-N_ν N^ν) for timelike vector + + if (not cmp::AlmostZero(norm)) { + Vector(i1, i2, i3, c_u0) = N_cntrv[0] / norm; + Vector(i1, i2, i3, c_u1) = N_cntrv[1] / norm; + Vector(i1, i2, i3, c_u2) = N_cntrv[2] / norm; + Vector(i1, i2, i3, c_u3) = N_cntrv[3] / norm; + } + } else { + raise::KernelError( + HERE, + "3D implementation of Normalize4VelocityByNorm_kernel called for non-3D"); + } + } + }; + + template + requires metric::traits::HasD && metric::traits::HasTransform + class Transform4VelocitySpatialToPhysical_kernel { + // Transforms spatial components of 4-velocity from coordinate to physical basis + // u^0 (Gamma/alpha) remains unchanged as it's unitless + ndfield_t Vector; + const unsigned short c_u1, c_u2, c_u3; + const M metric; + + public: + Transform4VelocitySpatialToPhysical_kernel(ndfield_t& vector, + unsigned short cu1, + unsigned short cu2, + unsigned short cu3, + const M& metric) + : Vector { vector } + , c_u1 { cu1 } + , c_u2 { cu2 } + , c_u3 { cu3 } + , metric { metric } { + raise::ErrorIf(c_u1 >= N or c_u2 >= N or c_u3 >= N, + "Invalid component index", + HERE); + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1) + HALF; + x_Code[1] = static_cast(i2) + HALF; + + vec_t u_cntrv { Vector(i1, i2, c_u1), + Vector(i1, i2, c_u2), + Vector(i1, i2, c_u3) }; + + vec_t u_phys { ZERO }; + // Transform spatial components from coordinate contravariant to physical contravariant + metric.template transform(x_Code, u_cntrv, u_phys); + + Vector(i1, i2, c_u1) = u_phys[0]; + Vector(i1, i2, c_u2) = u_phys[1]; + Vector(i1, i2, c_u3) = u_phys[2]; + } else { + raise::KernelError( + HERE, + "2D implementation of Transform4VelocitySpatialToPhysical_kernel called for non-2D"); + } + } + + Inline void operator()(index_t i1, index_t i2, index_t i3) const { + if constexpr (D == Dim::_3D) { + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1) + HALF; + x_Code[1] = static_cast(i2) + HALF; + x_Code[2] = static_cast(i3) + HALF; + + vec_t u_cntrv { Vector(i1, i2, i3, c_u1), + Vector(i1, i2, i3, c_u2), + Vector(i1, i2, i3, c_u3) }; + + vec_t u_phys { ZERO }; + // Transform spatial components from coordinate contravariant to physical contravariant + metric.template transform(x_Code, u_cntrv, u_phys); + + Vector(i1, i2, i3, c_u1) = u_phys[0]; + Vector(i1, i2, i3, c_u2) = u_phys[1]; + Vector(i1, i2, i3, c_u3) = u_phys[2]; + } else { + raise::KernelError( + HERE, + "3D implementation of Transform4VelocitySpatialToPhysical_kernel called for non-3D"); + } + } + }; + } // namespace kernel #endif // KERNELS_PARTICLE_MOMENTS_HPP diff --git a/src/output/fields.cpp b/src/output/fields.cpp index e6b86296f..c8f37c579 100644 --- a/src/output/fields.cpp +++ b/src/output/fields.cpp @@ -36,27 +36,44 @@ namespace out { raise::ErrorIf(id() == FldsID::A and S != SimEngine::GRPIC, "Output of A_phi not supported for non-GRPIC", HERE); - raise::ErrorIf(id() == FldsID::V and S == SimEngine::GRPIC, - "Output of bulk 3-vel not supported for GRPIC", - HERE); // determine the species and components to output if (is_moment()) { species = InterpretSpecies(name); } else { species = {}; } - if (is_field() || is_current() || id() == FldsID::V) { - // always write all the field/current/bulk vel components + if (is_field() || is_current()) { + // always write all the field/current components (3D) comp = { { 1 }, { 2 }, { 3 } }; + } else if (id() == FldsID::V) { + // bulk velocity: 3D for SR, 4D for GR + if (S == SimEngine::SRPIC) { + // SR: always 3-velocity components + comp = { { 1 }, { 2 }, { 3 } }; + } else { + // GR: 4-velocity components (Eckart frame) + // Check if user specified components + if (name.length() > 1) { + comp = InterpretComponents({ name.substr(1, 1) }); + // Validate component indices are in [0,3] for 4-vector + for (auto& c : comp) { + for (auto& idx : c) { + if (idx > 3) { + raise::Error("Invalid component index for 4-velocity (must be 0-3)", HERE); + } + } + } + } else { + // No component specified, output all 4 + comp = { { 0 }, { 1 }, { 2 }, { 3 } }; + } + } } else if (id() == FldsID::A) { // only write A3 comp = { { 3 } }; } else if (id() == FldsID::T) { // energy-momentum tensor comp = InterpretComponents({ name.substr(1, 1), name.substr(2, 1) }); - } else if (id() == FldsID::V) { - // energy-momentum tensor - comp = InterpretComponents({ name.substr(1, 1) }); } else { // scalar (Rho, divE, Custom, etc.) comp = {}; From 365b33d41b76f331a34624eff465e89abb01af18 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 4 Mar 2026 14:14:23 -0500 Subject: [PATCH 14/22] Added safety fallback to rest frame Eckart velocity for empty cells --- src/kernels/particle_moments.hpp | 69 +++++++++++++++++++++++++++----- 1 file changed, 58 insertions(+), 11 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 1dc59d6af..e4e9f980f 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -467,26 +467,50 @@ namespace kernel { coord_t x_Code { ZERO }; x_Code[0] = static_cast(i1) + HALF; x_Code[1] = static_cast(i2) + HALF; - + vec_t N_cntrv { ZERO }; N_cntrv[0] = Flux(i1, i2, c_u0); N_cntrv[1] = Flux(i1, i2, c_u1); N_cntrv[2] = Flux(i1, i2, c_u2); N_cntrv[3] = Flux(i1, i2, c_u3); - + + // rest frame for empty cells + if (cmp::AlmostZero(N_cntrv[0])) { + Vector(i1, i2, c_u0) = ONE; + Vector(i1, i2, c_u1) = ZERO; + Vector(i1, i2, c_u2) = ZERO; + Vector(i1, i2, c_u3) = ZERO; + return; + } + vec_t N_cov { ZERO }; // Compute N_i = g_ij N^j metric.template transform_4d(x_Code, N_cntrv, N_cov); - - // Compute N_ν N^ν = -(N^0)^2 + g_ij N^i N^j + + // Compute N_ν N^ν = g_μν N^μ N^ν (should be negative for timelike) real_t N_norm_sq { N_cov[0] * N_cntrv[0] + N_cov[1] * N_cntrv[1] + N_cov[2] * N_cntrv[2] + N_cov[3] * N_cntrv[3] }; + + // Set to rest frame for cells with insufficient particles + if (cmp::AlmostZero(N_norm_sq) || N_norm_sq >= ZERO) { + Vector(i1, i2, c_u0) = ONE; + Vector(i1, i2, c_u1) = ZERO; + Vector(i1, i2, c_u2) = ZERO; + Vector(i1, i2, c_u3) = ZERO; + return; + } + real_t norm = math::sqrt(-N_norm_sq); - + if (not cmp::AlmostZero(norm)) { Vector(i1, i2, c_u0) = N_cntrv[0] / norm; Vector(i1, i2, c_u1) = N_cntrv[1] / norm; Vector(i1, i2, c_u2) = N_cntrv[2] / norm; Vector(i1, i2, c_u3) = N_cntrv[3] / norm; + } else { + Vector(i1, i2, c_u0) = ONE; + Vector(i1, i2, c_u1) = ZERO; + Vector(i1, i2, c_u2) = ZERO; + Vector(i1, i2, c_u3) = ZERO; } } else { raise::KernelError( @@ -501,27 +525,50 @@ namespace kernel { x_Code[0] = static_cast(i1) + HALF; x_Code[1] = static_cast(i2) + HALF; x_Code[2] = static_cast(i3) + HALF; - + vec_t N_cntrv { ZERO }; N_cntrv[0] = Flux(i1, i2, i3, c_u0); // N^0 N_cntrv[1] = Flux(i1, i2, i3, c_u1); // N^1 N_cntrv[2] = Flux(i1, i2, i3, c_u2); // N^2 N_cntrv[3] = Flux(i1, i2, i3, c_u3); // N^3 - + + // Rest frame for empty cells + if (cmp::AlmostZero(N_cntrv[0])) { + Vector(i1, i2, i3, c_u0) = ONE; + Vector(i1, i2, i3, c_u1) = ZERO; + Vector(i1, i2, i3, c_u2) = ZERO; + Vector(i1, i2, i3, c_u3) = ZERO; + return; + } + vec_t N_cov { ZERO }; // Compute N_i = g_ij N^j metric.template transform_4d(x_Code, N_cntrv, N_cov); - - // Compute N_ν N^ν = -(N^0)^2 + g_ij N^i N^j + + // Compute N_ν N^ν = g_μν N^μ N^ν (should be negative for timelike) real_t N_norm_sq { N_cov[0] * N_cntrv[0] + N_cov[1] * N_cntrv[1] + N_cov[2] * N_cntrv[2] + N_cov[3] * N_cntrv[3] }; - + + // Set to rest frame for cells with insufficient particles + if (cmp::AlmostZero(N_norm_sq) || N_norm_sq >= ZERO) { + Vector(i1, i2, i3, c_u0) = ONE; + Vector(i1, i2, i3, c_u1) = ZERO; + Vector(i1, i2, i3, c_u2) = ZERO; + Vector(i1, i2, i3, c_u3) = ZERO; + return; + } + real_t norm = math::sqrt(-N_norm_sq); // sqrt(-N_ν N^ν) for timelike vector - + if (not cmp::AlmostZero(norm)) { Vector(i1, i2, i3, c_u0) = N_cntrv[0] / norm; Vector(i1, i2, i3, c_u1) = N_cntrv[1] / norm; Vector(i1, i2, i3, c_u2) = N_cntrv[2] / norm; Vector(i1, i2, i3, c_u3) = N_cntrv[3] / norm; + } else { + Vector(i1, i2, i3, c_u0) = ONE; + Vector(i1, i2, i3, c_u1) = ZERO; + Vector(i1, i2, i3, c_u2) = ZERO; + Vector(i1, i2, i3, c_u3) = ZERO; } } else { raise::KernelError( From f71fbbf7443eef4e2e4e501d1c41c92bf14ff4d1 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:56:08 -0400 Subject: [PATCH 15/22] removed unnec calculation --- src/kernels/particle_moments.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index e4e9f980f..54f181104 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -257,10 +257,6 @@ namespace kernel { u_cntrv_4d); // in GR: u^0 = Gamma/alpha u0 = u_cntrv_4d[0]; - metric.template transform( - x_Code, - { u_cntrv_4d[1], u_cntrv_4d[2], u_cntrv_4d[3] }, - u_Phys); // Deposit flux N^μ = mass * u^μ / u^0 if (c1 == 0) { // u^0 component From d55f170d2d277fb9bd04714d4934a62656b048e7 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:56:50 -0400 Subject: [PATCH 16/22] fixed missing NGHOSTS --- src/kernels/particle_moments.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 54f181104..e38d59c7e 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -461,8 +461,8 @@ namespace kernel { Inline void operator()(index_t i1, index_t i2) const { if constexpr (D == Dim::_2D) { coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1) + HALF; - x_Code[1] = static_cast(i2) + HALF; + x_Code[0] = COORD(i1) + HALF; + x_Code[1] = COORD(i2) + HALF; vec_t N_cntrv { ZERO }; N_cntrv[0] = Flux(i1, i2, c_u0); @@ -518,9 +518,9 @@ namespace kernel { Inline void operator()(index_t i1, index_t i2, index_t i3) const { if constexpr (D == Dim::_3D) { coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1) + HALF; - x_Code[1] = static_cast(i2) + HALF; - x_Code[2] = static_cast(i3) + HALF; + x_Code[0] = COORD(i1) + HALF; + x_Code[1] = COORD(i2) + HALF; + x_Code[2] = COORD(i3) + HALF; vec_t N_cntrv { ZERO }; N_cntrv[0] = Flux(i1, i2, i3, c_u0); // N^0 @@ -602,8 +602,8 @@ namespace kernel { Inline void operator()(index_t i1, index_t i2) const { if constexpr (D == Dim::_2D) { coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1) + HALF; - x_Code[1] = static_cast(i2) + HALF; + x_Code[0] = COORD(i1) + HALF; + x_Code[1] = COORD(i2) + HALF; vec_t u_cntrv { Vector(i1, i2, c_u1), Vector(i1, i2, c_u2), @@ -626,9 +626,9 @@ namespace kernel { Inline void operator()(index_t i1, index_t i2, index_t i3) const { if constexpr (D == Dim::_3D) { coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1) + HALF; - x_Code[1] = static_cast(i2) + HALF; - x_Code[2] = static_cast(i3) + HALF; + x_Code[0] = COORD(i1) + HALF; + x_Code[1] = COORD(i2) + HALF; + x_Code[2] = COORD(i3) + HALF; vec_t u_cntrv { Vector(i1, i2, i3, c_u1), Vector(i1, i2, i3, c_u2), From 6e54a61c465a95ca7e23ed906cf3025e06a87137 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:57:32 -0400 Subject: [PATCH 17/22] fallback to ZAMO when no enough particles --- src/kernels/particle_moments.hpp | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index e38d59c7e..e5f2fbfeb 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -470,12 +470,19 @@ namespace kernel { N_cntrv[2] = Flux(i1, i2, c_u2); N_cntrv[3] = Flux(i1, i2, c_u3); - // rest frame for empty cells - if (cmp::AlmostZero(N_cntrv[0])) { - Vector(i1, i2, c_u0) = ONE; - Vector(i1, i2, c_u1) = ZERO; + // ZAMO fallback for empty or pathological cells + const auto zamo_fallback_2d = [&]() { + const real_t al = metric.alpha(x_Code); + const real_t b1 = metric.beta1(x_Code); + Vector(i1, i2, c_u0) = ONE / al; + Vector(i1, i2, c_u1) = -b1 / al; Vector(i1, i2, c_u2) = ZERO; Vector(i1, i2, c_u3) = ZERO; + }; + + // rest frame for empty cells + if (cmp::AlmostZero(N_cntrv[0])) { + zamo_fallback_2d(); return; } @@ -486,12 +493,9 @@ namespace kernel { // Compute N_ν N^ν = g_μν N^μ N^ν (should be negative for timelike) real_t N_norm_sq { N_cov[0] * N_cntrv[0] + N_cov[1] * N_cntrv[1] + N_cov[2] * N_cntrv[2] + N_cov[3] * N_cntrv[3] }; - // Set to rest frame for cells with insufficient particles + // ZAMO fallback for cells with insufficient particles if (cmp::AlmostZero(N_norm_sq) || N_norm_sq >= ZERO) { - Vector(i1, i2, c_u0) = ONE; - Vector(i1, i2, c_u1) = ZERO; - Vector(i1, i2, c_u2) = ZERO; - Vector(i1, i2, c_u3) = ZERO; + zamo_fallback_2d(); return; } From 8974d069bb442694e2aa78de1452cc8602e6dd2e Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Thu, 12 Mar 2026 10:51:31 -0400 Subject: [PATCH 18/22] removed unused variable --- src/kernels/particle_moments.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index e5f2fbfeb..4bc447d84 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -237,7 +237,6 @@ namespace kernel { } else { // GR: Eckart frame flux N^μ = m * u^μ / u^0 static_assert(D != Dim::_1D, "GRPIC 1D"); - vec_t u_Phys { ZERO }; real_t u0 { ZERO }; coord_t x_Code { ZERO }; x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); From 331d641a4ccca94a8452c3200efced4059f66bb3 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Thu, 12 Mar 2026 15:44:05 -0400 Subject: [PATCH 19/22] simplified Tmunu and bulk velocity calculation --- src/kernels/particle_moments.hpp | 249 ++++++++++++++++--------------- 1 file changed, 126 insertions(+), 123 deletions(-) diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 4bc447d84..5b98ecb45 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -131,139 +131,142 @@ namespace kernel { } } - Inline void operator()(index_t p) const { - if (tag(p) == ParticleTag::dead) { - return; - } - real_t coeff { ZERO }; - if constexpr (F == FldsID::T) { - real_t u0 { ZERO }; - // for stress-energy tensor - vec_t u_Phys { ZERO }; - if constexpr (S == SimEngine::SRPIC) { - // SR - // stress-energy tensor for SR is computed in the tetrad (hatted) basis - if constexpr (M::CoordType == Coord::Cart) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); - } else { - static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D"); - coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); - if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); - } else { - x_Code[2] = phi(p); - } - metric.template transform_xyz( - x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); - } - if (mass == ZERO) { - u0 = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); - } else { - u0 = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); - } + Inline auto computeStressEnergyComponent(index_t p) const -> real_t { + real_t u0 { ZERO }; + vec_t u_Phys { ZERO }; + if constexpr (S == SimEngine::SRPIC) { + // stress-energy tensor for SR is computed in the tetrad (hatted) basis + if constexpr (M::CoordType == Coord::Cart) { + u_Phys[0] = ux1(p); + u_Phys[1] = ux2(p); + u_Phys[2] = ux3(p); } else { - // GR - // stress-energy tensor for GR is computed in contravariant basis - static_assert(D != Dim::_1D, "GRPIC 1D"); - coord_t x_Code { ZERO }; + static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D"); + coord_t x_Code { ZERO }; x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); if constexpr (D == Dim::_3D) { x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + } else { + x_Code[2] = phi(p); } - // raise full covariant 4-vector to get correct contravariant u^i - // u^i != h^{ij} u_j - const real_t u_0_cov { metric.u_0(x_Code, - { ux1(p), ux2(p), ux3(p) }, - (mass == ZERO) ? ZERO : ONE) }; - vec_t u_cntrv_4d { ZERO }; - metric.template transform_4d( - x_Code, - { u_0_cov, ux1(p), ux2(p), ux3(p) }, - u_cntrv_4d); - // in GR: u^0 = Gamma/alpha - u0 = u_cntrv_4d[0]; - metric.template transform( - x_Code, - { u_cntrv_4d[1], u_cntrv_4d[2], u_cntrv_4d[3] }, - u_Phys); + metric.template transform_xyz(x_Code, + { ux1(p), ux2(p), ux3(p) }, + u_Phys); } - // compute the corresponding moment - // T^μν = m * u^0 * v^{c1} * v^{c2}, - // where v^0 = 1, v^i = u^i / u^0 (coordinate 3-velocity) - coeff = ((mass == ZERO) ? ONE : mass) * u0; -#pragma unroll - for (const auto& c : { c1, c2 }) { - if (c != 0) { - coeff *= u_Phys[c - 1] / u0; - } + u0 = (mass == ZERO) + ? (NORM(u_Phys[0], u_Phys[1], u_Phys[2])) + : (math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]))); + } else if constexpr (S == SimEngine::GRPIC) { + // stress-energy tensor for GR is computed in contravariant basis + static_assert(D != Dim::_1D, "GRPIC 1D"); + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); + x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + if constexpr (D == Dim::_3D) { + x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); } + // raise full covariant 4-vector to get correct contravariant u^i + // u^i != h^{ij} u_j + const real_t u_0_cov { metric.u_0(x_Code, + { ux1(p), ux2(p), ux3(p) }, + (mass == ZERO) ? ZERO : ONE) }; + vec_t u_cntrv_4d { ZERO }; + metric.template transform_4d( + x_Code, + { u_0_cov, ux1(p), ux2(p), ux3(p) }, + u_cntrv_4d); + // in GR: u^0 = Gamma/alpha + u0 = u_cntrv_4d[0]; + metric.template transform( + x_Code, + { u_cntrv_4d[1], u_cntrv_4d[2], u_cntrv_4d[3] }, + u_Phys); + } else { + raise::KernelError( + HERE, + "computeStressEnergyComponent called for non-SRPIC/GRPIC"); + } + auto T_component = (mass == ZERO ? ONE : mass) / u0; + for (const auto& c : { c1, c2 }) { + if (c > 0) { + T_component *= u_Phys[c - 1]; + } else { + T_component *= u0; + } + } + return T_component; + } + + Inline auto computeBulk3VelocityTimesMass(index_t p) const -> real_t { + real_t u0 { ZERO }; + // for bulk 3vel (tetrad basis) + vec_t u_Phys { ZERO }; + if constexpr (M::CoordType == Coord::Cart) { + u_Phys[0] = ux1(p); + u_Phys[1] = ux2(p); + u_Phys[2] = ux3(p); + } else { + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); + x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + if constexpr (D == Dim::_3D) { + x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + } else { + x_Code[2] = phi(p); + } + metric.template transform_xyz(x_Code, + { ux1(p), ux2(p), ux3(p) }, + u_Phys); + } + if (mass == ZERO) { + u0 = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); + } else { + u0 = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); + } + return (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / u0; + } + + Inline auto computeEckartVelocityFluxComponent(index_t p) const -> real_t { + // GR: Eckart frame flux N^μ = m * u^μ / u^0 + static_assert(D != Dim::_1D, "GRPIC 1D"); + coord_t x_Code { ZERO }; + x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); + x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + if constexpr (D == Dim::_3D) { + x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + } + // raise full covariant 4-vector to get correct contravariant u^μ + // u^i != h^{ij} u_j + const real_t u_0_cov { metric.u_0(x_Code, + { ux1(p), ux2(p), ux3(p) }, + (mass == ZERO) ? ZERO : ONE) }; + vec_t u_cntrv_4d { ZERO }; + metric.template transform_4d( + x_Code, + { u_0_cov, ux1(p), ux2(p), ux3(p) }, + u_cntrv_4d); + const real_t u0 { u_cntrv_4d[0] }; + // Deposit flux N^μ = mass * u^μ / u^0 + if (c1 == 0) { + return (mass == ZERO ? ONE : mass); + } else { + return (mass == ZERO ? ONE : mass) * u_cntrv_4d[c1] / u0; + } + } + + Inline void operator()(index_t p) const { + if (tag(p) == ParticleTag::dead) { + return; + } + real_t coeff { ZERO }; + if constexpr (F == FldsID::T) { + coeff = computeStressEnergyComponent(p); } else if constexpr (F == FldsID::V) { - if constexpr (S == SimEngine::SRPIC) { - // SR: bulk velocity - real_t gamma { ZERO }; - // for bulk 3vel (tetrad basis) - vec_t u_Phys { ZERO }; - if constexpr (M::CoordType == Coord::Cart) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); - } else { - coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); - if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); - } else { - x_Code[2] = phi(p); - } - metric.template transform_xyz(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); - } - if (mass == ZERO) { - gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); - } else { - gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); - } - // compute the corresponding moment - coeff = (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / gamma; + if constexpr (S == SimEngine::GRPIC) { + coeff = computeEckartVelocityFluxComponent(p); } else { - // GR: Eckart frame flux N^μ = m * u^μ / u^0 - static_assert(D != Dim::_1D, "GRPIC 1D"); - real_t u0 { ZERO }; - coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); - if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); - } - // raise full covariant 4-vector to get correct contravariant u^i - // u^i != h^{ij} u_j - const real_t u_0_cov { metric.u_0(x_Code, - { ux1(p), ux2(p), ux3(p) }, - (mass == ZERO) ? ZERO : ONE) }; - vec_t u_cntrv_4d { ZERO }; - metric.template transform_4d( - x_Code, - { u_0_cov, ux1(p), ux2(p), ux3(p) }, - u_cntrv_4d); - // in GR: u^0 = Gamma/alpha - u0 = u_cntrv_4d[0]; - // Deposit flux N^μ = mass * u^μ / u^0 - if (c1 == 0) { - // u^0 component - coeff = (mass == ZERO ? ONE : mass); - } else { - // u^i component: mass * u^i / u^0 - coeff = (mass == ZERO ? ONE : mass) * u_cntrv_4d[c1] / u0; - } + coeff = computeBulk3VelocityTimesMass(p); } } else { // for other cases, use the `contrib` defined above From 342c8008e80190799fca9ea6ad968c5914fdb4c9 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Thu, 12 Mar 2026 16:12:28 -0400 Subject: [PATCH 20/22] always output 4 components of bulk vel in GR --- src/output/fields.cpp | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/output/fields.cpp b/src/output/fields.cpp index c8f37c579..7e3d05cab 100644 --- a/src/output/fields.cpp +++ b/src/output/fields.cpp @@ -51,22 +51,8 @@ namespace out { // SR: always 3-velocity components comp = { { 1 }, { 2 }, { 3 } }; } else { - // GR: 4-velocity components (Eckart frame) - // Check if user specified components - if (name.length() > 1) { - comp = InterpretComponents({ name.substr(1, 1) }); - // Validate component indices are in [0,3] for 4-vector - for (auto& c : comp) { - for (auto& idx : c) { - if (idx > 3) { - raise::Error("Invalid component index for 4-velocity (must be 0-3)", HERE); - } - } - } - } else { - // No component specified, output all 4 - comp = { { 0 }, { 1 }, { 2 }, { 3 } }; - } + // GR: always output all 4 Eckart velocity components (u^0, u^1, u^2, u^3) + comp = { { 0 }, { 1 }, { 2 }, { 3 } }; } } else if (id() == FldsID::A) { // only write A3 From 36afb1da0384a3781d56913d5a846728a5abcd12 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Thu, 12 Mar 2026 16:19:19 -0400 Subject: [PATCH 21/22] removed useless blocks since "V" always outputs 3 components for SR and 4 components for GR --- src/framework/domain/output.cpp | 44 ++++++--------------------------- 1 file changed, 7 insertions(+), 37 deletions(-) diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 9601b1ddb..dadd7a514 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -486,25 +486,6 @@ namespace ntt { {}, local_domain->fields.bckp, c); - } else if (fld.id() == FldsID::V) { - if constexpr (S == SimEngine::SRPIC) { - ComputeMoments(params, - local_domain->mesh, - local_domain->species, - fld.species, - fld.comp[0], - local_domain->fields.bckp, - c); - } else if constexpr (S == SimEngine::GRPIC) { - // Single component of 4-velocity (user requested e.g., "v_0") - ComputeMoments(params, - local_domain->mesh, - local_domain->species, - fld.species, - fld.comp[0], - local_domain->fields.bckp, - c); - } } else { raise::Error("Wrong moment requested for output", HERE); } @@ -576,24 +557,13 @@ namespace ntt { raise::ErrorIf(fld.comp[i].size() != 1, "Wrong # of components requested for 3vel", HERE); - if constexpr (S == SimEngine::SRPIC) { - ComputeMoments(params, - local_domain->mesh, - local_domain->species, - fld.species, - fld.comp[i], - local_domain->fields.bckp, - c); - } else if constexpr (S == SimEngine::GRPIC) { - // Single component from 3-vector spec (legacy), deposit as part of 4-vector - ComputeMoments(params, - local_domain->mesh, - local_domain->species, - fld.species, - fld.comp[i], - local_domain->fields.bckp, - c); - } + ComputeMoments(params, + local_domain->mesh, + local_domain->species, + fld.species, + fld.comp[i], + local_domain->fields.bckp, + c); } else { raise::Error("Wrong moment requested for output", HERE); } From cea4557e8e06107f8bb2db16fa64fb74248d0f34 Mon Sep 17 00:00:00 2001 From: StaticObserver Date: Fri, 20 Mar 2026 18:43:27 +0800 Subject: [PATCH 22/22] metric tests for 4d related functions --- src/metrics/tests/CMakeLists.txt | 1 + src/metrics/tests/ks-qks-4d.cpp | 542 +++++++++++++++++++++++++++++++ 2 files changed, 543 insertions(+) create mode 100644 src/metrics/tests/ks-qks-4d.cpp diff --git a/src/metrics/tests/CMakeLists.txt b/src/metrics/tests/CMakeLists.txt index 0d661c318..f1c43ccb7 100644 --- a/src/metrics/tests/CMakeLists.txt +++ b/src/metrics/tests/CMakeLists.txt @@ -29,3 +29,4 @@ gen_test(coord_trans) gen_test(sph-qsph) gen_test(ks-qks) gen_test(sr-cart-sph) +gen_test(ks-qks-4d) diff --git a/src/metrics/tests/ks-qks-4d.cpp b/src/metrics/tests/ks-qks-4d.cpp new file mode 100644 index 000000000..d46a9b676 --- /dev/null +++ b/src/metrics/tests/ks-qks-4d.cpp @@ -0,0 +1,542 @@ +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/comparators.h" + +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/qkerr_schild.h" + +#include +#include +#include +#include + +void errorIf(bool condition, const std::string& message) { + if (condition) { + throw std::runtime_error(message); + } +} + +inline static constexpr auto epsilon = std::numeric_limits::epsilon(); + +template +Inline auto equal(const vec_t& a, + const vec_t& b, + const char* msg, + real_t acc = ONE) -> bool { + const auto eps = epsilon * acc; + for (auto d { 0u }; d < D; ++d) { + if (not cmp::AlmostEqual(a[d], b[d], eps)) { + Kokkos::printf("%s [%d]: %.12e != %.12e\n", msg, d, a[d], b[d]); + return false; + } + } + return true; +} + +Inline auto almostZero(real_t val, real_t scale, real_t acc) -> bool { + return math::abs(val) <= scale * epsilon * acc; +} + +template +Inline void unravel(std::size_t idx, + tuple_t& ijk, + const tuple_t& res) { + for (auto d { 0u }; d < D; ++d) { + ijk[d] = idx % res[d]; + idx /= res[d]; + } +} + +/** + * Test g^{mu alpha} * g_{alpha nu} = delta^mu_nu + */ +template +void testGInverse(const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params) { + static_assert(M::Dim == 2, "Dim != 2"); + + M metric(res, ext, params); + tuple_t res_tup; + std::size_t npts = 1; + for (auto d = 0; d < M::Dim; ++d) { + res_tup[d] = res[d]; + npts *= res[d]; + } + + Kokkos::parallel_for( + "g_inverse", + npts, + Lambda(index_t n) { + tuple_t idx; + unravel(n, idx, res_tup); + coord_t x { ZERO }; + for (auto d { 0u }; d < M::Dim; ++d) { + x[d] = (real_t)(idx[d]) + HALF; + } + + const auto g00 = metric.template g<0, 0>(x); + const auto g01 = metric.template g<0, 1>(x); + const auto g11 = metric.template g<1, 1>(x); + const auto g22 = metric.template g<2, 2>(x); + const auto g33 = metric.template g<3, 3>(x); + const auto g13 = metric.template g<1, 3>(x); + + const auto gl00 = metric.template g_<0, 0>(x); + const auto gl01 = metric.template g_<0, 1>(x); + const auto gl03 = metric.template g_<0, 3>(x); + const auto gl11 = metric.template g_<1, 1>(x); + const auto gl22 = metric.template g_<2, 2>(x); + const auto gl33 = metric.template g_<3, 3>(x); + const auto gl13 = metric.template g_<1, 3>(x); + + // g^{0 alpha} g_{alpha 0} = 1 + const auto d00 = g00 * gl00 + g01 * gl01; + // g^{1 alpha} g_{alpha 1} = 1 + const auto d11 = g01 * gl01 + g11 * gl11 + g13 * gl13; + // g^{2 alpha} g_{alpha 2} = 1 + const auto d22 = g22 * gl22; + // g^{3 alpha} g_{alpha 3} = 1 + const auto d33 = g13 * gl13 + g33 * gl33; + + vec_t diag { d00, d11, d22, d33 }; + vec_t diag_expect { ONE, ONE, ONE, ONE }; + if (not equal(diag, diag_expect, "g_inv_diag", acc)) { + Kokkos::abort("g inverse: diagonal != identity"); + } + + // off-diagonal should be zero + // g^{0 alpha} g_{alpha 1} = 0 + const auto od01 = g00 * gl01 + g01 * gl11; + // g^{0 alpha} g_{alpha 3} = 0 + const auto od03 = g00 * gl03 + g01 * gl13; + // g^{1 alpha} g_{alpha 0} = 0 + const auto od10 = g01 * gl00 + g11 * gl01 + g13 * gl03; + // g^{1 alpha} g_{alpha 3} = 0 + const auto od13 = g01 * gl03 + g11 * gl13 + g13 * gl33; + // g^{3 alpha} g_{alpha 0} = 0 + const auto od30 = g13 * gl01 + g33 * gl03; + // g^{3 alpha} g_{alpha 1} = 0 + const auto od31 = g13 * gl11 + g33 * gl13; + + const auto scale = math::max(math::abs(gl11), math::abs(gl33)); + if (not almostZero(od01, scale, acc) || + not almostZero(od03, scale, acc) || + not almostZero(od10, scale, acc) || + not almostZero(od13, scale, acc) || + not almostZero(od30, scale, acc) || + not almostZero(od31, scale, acc)) { + Kokkos::printf("off-diag: %e %e %e %e %e %e\n", + od01, od03, od10, od13, od30, od31); + Kokkos::abort("g inverse: off-diagonal != 0"); + } + }); +} + +/** + * Test spatial components: g_{ij} == h_{ij} + * Valid for all KS-family metrics (KS0, KS, QKS). + */ +template +void testGSpatial(const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params) { + static_assert(M::Dim == 2, "Dim != 2"); + + M metric(res, ext, params); + tuple_t res_tup; + std::size_t npts = 1; + for (auto d = 0; d < M::Dim; ++d) { + res_tup[d] = res[d]; + npts *= res[d]; + } + + Kokkos::parallel_for( + "g_spatial", + npts, + Lambda(index_t n) { + tuple_t idx; + unravel(n, idx, res_tup); + coord_t x_Code { ZERO }; + for (auto d { 0u }; d < M::Dim; ++d) { + x_Code[d] = (real_t)(idx[d]) + HALF; + } + + // g_{ij} == h_{ij} for spatial components + vec_t g_ij { + metric.template g_<1, 1>(x_Code), + metric.template g_<2, 2>(x_Code), + metric.template g_<3, 3>(x_Code), + metric.template g_<1, 3>(x_Code) + }; + vec_t h_ij_expect { + metric.template h_<1, 1>(x_Code), + metric.template h_<2, 2>(x_Code), + metric.template h_<3, 3>(x_Code), + metric.template h_<1, 3>(x_Code) + }; + if (not equal(g_ij, h_ij_expect, "g_ij==h_ij", acc)) { + Kokkos::abort("g_spatial: g_{ij} != h_{ij}"); + } + }); +} + +/** + * Test g_00 and g^00 against analytical Kerr formulas: z = 2r/Sigma. + * Only valid for KerrSchild and QKerrSchild (NOT KerrSchild0). + */ +template +void testG00Analytical(const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params) { + static_assert(M::Dim == 2, "Dim != 2"); + + M metric(res, ext, params); + tuple_t res_tup; + std::size_t npts = 1; + for (auto d = 0; d < M::Dim; ++d) { + res_tup[d] = res[d]; + npts *= res[d]; + } + + const auto a_spin = metric.spin(); + + Kokkos::parallel_for( + "g00_analytical", + npts, + Lambda(index_t n) { + tuple_t idx; + unravel(n, idx, res_tup); + coord_t x_Code { ZERO }; + coord_t x_Phys { ZERO }; + for (auto d { 0u }; d < M::Dim; ++d) { + x_Code[d] = (real_t)(idx[d]) + HALF; + } + + metric.template convert(x_Code, x_Phys); + const auto r = x_Phys[0]; + const auto theta = x_Phys[1]; + const auto Sigma = SQR(r) + SQR(a_spin * math::cos(theta)); + const auto z_val = TWO * r / Sigma; + + vec_t g00_code { metric.template g_<0, 0>(x_Code) }; + vec_t g00_expect { -(ONE - z_val) }; + if (not equal(g00_code, g00_expect, "g_00", acc)) { + Kokkos::abort("g00_analytical: g_00"); + } + + vec_t ginv00_code { metric.template g<0, 0>(x_Code) }; + vec_t ginv00_expect { -(ONE + z_val) }; + if (not equal(ginv00_code, ginv00_expect, "g^00", acc)) { + Kokkos::abort("g00_analytical: g^00"); + } + }); +} + +/** + * Test u_0: verify g^{mu nu} u_mu u_nu + norm = 0 + * The u_0 function solves: g^{mu nu} u_mu u_nu = -norm + * For massive particles (g^{mu nu} u_mu u_nu = -1): norm = 1 + * For null (g^{mu nu} u_mu u_nu = 0): norm = 0 + */ +template +void testU0(const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params) { + static_assert(M::Dim == 2, "Dim != 2"); + + M metric(res, ext, params); + tuple_t res_tup; + std::size_t npts = 1; + for (auto d = 0; d < M::Dim; ++d) { + res_tup[d] = res[d]; + npts *= res[d]; + } + + Kokkos::parallel_for( + "u0_normalization", + npts, + Lambda(index_t n) { + tuple_t idx; + unravel(n, idx, res_tup); + coord_t x { ZERO }; + for (auto d { 0u }; d < M::Dim; ++d) { + x[d] = (real_t)(idx[d]) + HALF; + } + + const real_t norm = ONE; + + // Test with zero spatial velocity (covariant) + { + vec_t u_i { ZERO, ZERO, ZERO }; + const auto u0 = metric.u_0(x, u_i, norm); + + const auto g00 = metric.template g<0, 0>(x); + const auto residual = g00 * SQR(u0) + norm; + if (not almostZero(residual, norm, acc)) { + Kokkos::printf("u0_zero: g^00*u0^2+norm = %.12e (u0=%.12e)\n", + residual, u0); + Kokkos::abort("u_0: zero velocity normalization failed"); + } + } + + // Test with non-zero spatial velocity + { + vec_t u_i { (real_t)0.1, (real_t)0.05, (real_t)0.02 }; + const auto u0 = metric.u_0(x, u_i, norm); + + const auto g00 = metric.template g<0, 0>(x); + const auto g01 = metric.template g<0, 1>(x); + const auto g11 = metric.template g<1, 1>(x); + const auto g22 = metric.template g<2, 2>(x); + const auto g33 = metric.template g<3, 3>(x); + const auto g13 = metric.template g<1, 3>(x); + + const auto contraction = g00 * SQR(u0) + + TWO * g01 * u0 * u_i[0] + + g11 * SQR(u_i[0]) + + g22 * SQR(u_i[1]) + + g33 * SQR(u_i[2]) + + TWO * g13 * u_i[0] * u_i[2]; + const auto residual = contraction + norm; + + if (not almostZero(residual, norm, acc)) { + Kokkos::printf("u0_nonzero: g^uv u_u u_v + norm = %.12e\n", residual); + Kokkos::abort("u_0: nonzero velocity normalization failed"); + } + } + + // Test with norm = 0 (null vector) + { + vec_t u_i { (real_t)0.3, ZERO, ZERO }; + const auto u0 = metric.u_0(x, u_i, ZERO); + + const auto g00 = metric.template g<0, 0>(x); + const auto g01 = metric.template g<0, 1>(x); + const auto g11 = metric.template g<1, 1>(x); + + const auto contraction = g00 * SQR(u0) + + TWO * g01 * u0 * u_i[0] + + g11 * SQR(u_i[0]); + + if (not almostZero(contraction, math::abs(g00), acc)) { + Kokkos::printf("u0_null: g^uv u_u u_v = %.12e, expect 0\n", + contraction); + Kokkos::abort("u_0: null normalization failed"); + } + } + }); +} + +/** + * Test transform_4d: D->U->D and U->D->U roundtrip identity + */ +template +void testTransform4d(const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params) { + static_assert(M::Dim == 2, "Dim != 2"); + + M metric(res, ext, params); + tuple_t res_tup; + std::size_t npts = 1; + for (auto d = 0; d < M::Dim; ++d) { + res_tup[d] = res[d]; + npts *= res[d]; + } + + Kokkos::parallel_for( + "transform_4d", + npts, + Lambda(index_t n) { + tuple_t idx; + unravel(n, idx, res_tup); + coord_t x { ZERO }; + for (auto d { 0u }; d < M::Dim; ++d) { + x[d] = (real_t)(idx[d]) + HALF; + } + + // D -> U -> D roundtrip + { + vec_t v_cov { (real_t)-1.5, (real_t)0.3, (real_t)0.1, (real_t)0.2 }; + vec_t v_cntrv { ZERO }; + vec_t v_cov_back { ZERO }; + + metric.template transform_4d(x, v_cov, v_cntrv); + metric.template transform_4d(x, v_cntrv, v_cov_back); + + if (not equal(v_cov_back, v_cov, "D->U->D", acc)) { + Kokkos::abort("transform_4d: D->U->D roundtrip failed"); + } + } + + // U -> D -> U roundtrip + { + vec_t v_cntrv { (real_t)1.0, (real_t)0.5, (real_t)-0.2, (real_t)0.1 }; + vec_t v_cov { ZERO }; + vec_t v_cntrv_back { ZERO }; + + metric.template transform_4d(x, v_cntrv, v_cov); + metric.template transform_4d(x, v_cov, v_cntrv_back); + + if (not equal(v_cntrv_back, v_cntrv, "U->D->U", acc)) { + Kokkos::abort("transform_4d: U->D->U roundtrip failed"); + } + } + + // Verify transform_4d is consistent with direct g^{mu nu} multiplication + { + vec_t v_cov { (real_t)-2.0, (real_t)0.7, (real_t)0.3, (real_t)0.4 }; + vec_t v_cntrv_tr { ZERO }; + metric.template transform_4d(x, v_cov, v_cntrv_tr); + + const auto g00 = metric.template g<0, 0>(x); + const auto g01 = metric.template g<0, 1>(x); + const auto g11 = metric.template g<1, 1>(x); + const auto g22 = metric.template g<2, 2>(x); + const auto g33 = metric.template g<3, 3>(x); + const auto g13 = metric.template g<1, 3>(x); + + vec_t v_cntrv_direct { + g00 * v_cov[0] + g01 * v_cov[1], + g01 * v_cov[0] + g11 * v_cov[1] + g13 * v_cov[3], + g22 * v_cov[2], + g13 * v_cov[1] + g33 * v_cov[3] + }; + + if (not equal(v_cntrv_tr, v_cntrv_direct, "D->U consistency", acc)) { + Kokkos::abort("transform_4d: D->U inconsistent with g^{mu nu}"); + } + } + + // Verify transform_4d is consistent with direct g_{mu nu} multiplication + { + vec_t v_cntrv { (real_t)1.2, (real_t)-0.4, (real_t)0.6, (real_t)0.15 }; + vec_t v_cov_tr { ZERO }; + metric.template transform_4d(x, v_cntrv, v_cov_tr); + + const auto gl00 = metric.template g_<0, 0>(x); + const auto gl01 = metric.template g_<0, 1>(x); + const auto gl03 = metric.template g_<0, 3>(x); + const auto gl11 = metric.template g_<1, 1>(x); + const auto gl22 = metric.template g_<2, 2>(x); + const auto gl33 = metric.template g_<3, 3>(x); + const auto gl13 = metric.template g_<1, 3>(x); + + vec_t v_cov_direct { + gl00 * v_cntrv[0] + gl01 * v_cntrv[1] + gl03 * v_cntrv[3], + gl01 * v_cntrv[0] + gl11 * v_cntrv[1] + gl13 * v_cntrv[3], + gl22 * v_cntrv[2], + gl03 * v_cntrv[0] + gl13 * v_cntrv[1] + gl33 * v_cntrv[3] + }; + + if (not equal(v_cov_tr, v_cov_direct, "U->D consistency", acc)) { + Kokkos::abort("transform_4d: U->D inconsistent with g_{mu nu}"); + } + } + }); +} + +template +void runAllTests(const char* label, + const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params, + bool has_kerr_g00 = false) { + std::cout << "[" << label << "] testGInverse..." << std::endl; + testGInverse(res, ext, acc, params); + std::cout << " passed." << std::endl; + + std::cout << "[" << label << "] testGSpatial..." << std::endl; + testGSpatial(res, ext, acc, params); + std::cout << " passed." << std::endl; + + if (has_kerr_g00) { + std::cout << "[" << label << "] testG00Analytical..." << std::endl; + testG00Analytical(res, ext, acc, params); + std::cout << " passed." << std::endl; + } + + std::cout << "[" << label << "] testU0..." << std::endl; + testU0(res, ext, acc, params); + std::cout << " passed." << std::endl; + + std::cout << "[" << label << "] testTransform4d..." << std::endl; + testTransform4d(res, ext, acc, params); + std::cout << " passed." << std::endl; +} + +auto main(int argc, char* argv[]) -> int { + Kokkos::initialize(argc, argv); + + try { + using namespace ntt; + using namespace metric; + + const auto res = std::vector { 32, 42 }; + const auto ext = boundaries_t { + { 0.8, 10.0 }, + { 0.0, constant::PI } + }; + const real_t acc = static_cast(1e4); + + // ---- KerrSchild0 (a=0, flat spacetime in spherical coords) ---- + std::cout << "=== KerrSchild0 ===" << std::endl; + runAllTests>("KS0", res, ext, acc, {}); + + // ---- KerrSchild (full Kerr, spherical KS coords) ---- + std::cout << "=== KerrSchild (a=0.8) ===" << std::endl; + const auto ks_params_08 = std::map { + { "a", (real_t)0.8 } + }; + runAllTests>("KS a=0.8", res, ext, acc, + ks_params_08, true); + + std::cout << "=== KerrSchild (a=0.95) ===" << std::endl; + const auto ks_params_095 = std::map { + { "a", (real_t)0.95 } + }; + runAllTests>( + "KS a=0.95", { 64, 54 }, { { 0.8, 50.0 }, { 0.0, constant::PI } }, + acc, ks_params_095, true); + + // ---- QKerrSchild (h=0, quasi-spherical reduces to spherical) ---- + std::cout << "=== QKerrSchild (a=0.8, h=0) ===" << std::endl; + const auto qks_h0 = std::map { + { "r0", -TWO }, { "h", ZERO }, { "a", (real_t)0.8 } + }; + runAllTests>("QKS h=0", res, ext, acc, + qks_h0, true); + + // ---- QKerrSchild (h=0.25, with angular stretching) ---- + std::cout << "=== QKerrSchild (a=0.8, h=0.25) ===" << std::endl; + const auto qks_h025 = std::map { + { "r0", -TWO }, { "h", (real_t)0.25 }, { "a", (real_t)0.8 } + }; + runAllTests>("QKS h=0.25", res, ext, acc, + qks_h025, true); + + // ---- QKerrSchild (low spin) ---- + std::cout << "=== QKerrSchild (a=0.3, h=0.1) ===" << std::endl; + const auto qks_low = std::map { + { "r0", -TWO }, { "h", (real_t)0.1 }, { "a", (real_t)0.3 } + }; + runAllTests>("QKS a=0.3", res, ext, acc, + qks_low, true); + + } catch (std::exception& e) { + std::cerr << e.what() << std::endl; + Kokkos::finalize(); + return 1; + } + Kokkos::finalize(); + return 0; +}