diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 1d29c4102..dadd7a514 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -486,18 +486,6 @@ namespace ntt { {}, local_domain->fields.bckp, c); - } else if (fld.id() == FldsID::V) { - if constexpr (S != SimEngine::GRPIC) { - ComputeMoments(params, - local_domain->mesh, - local_domain->species, - fld.species, - 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); } @@ -569,17 +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 { - raise::Error("Bulk velocity not supported for GRPIC", HERE); - } + 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 +662,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 023eb55a5..5b98ecb45 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"); @@ -133,85 +131,17 @@ 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 energy { 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) { - energy = 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])); - } - } 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); - 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); - } - metric.template transform(x_Code, u_Cntrv, u_Phys); - } - // compute the corresponding moment - coeff = ONE / energy; -#pragma unroll - for (const auto& c : { c1, c2 }) { - if (c == 0) { - coeff *= energy; - } else { - coeff *= ((mass == ZERO) ? ONE : mass) * u_Phys[c - 1]; - } - } - } else if constexpr (F == FldsID::V) { - real_t gamma { ZERO }; - // for bulk 3vel (tetrad basis) - vec_t u_Phys { ZERO }; + 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 { + 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)); @@ -224,76 +154,124 @@ namespace kernel { { 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])); + 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)); } - // compute the corresponding moment - coeff = (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / gamma; + // 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 { - // for other cases, use the `contrib` defined above - coeff = contrib; + 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; + } - 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])); - } + 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 { - // 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); + x_Code[2] = phi(p); } - // compute the corresponding moment - coeff = u_Phys[c1 - 1] / gamma; + 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::GRPIC) { + coeff = computeEckartVelocityFluxComponent(p); + } else { + coeff = computeBulk3VelocityTimesMass(p); + } + } else { + // for other cases, use the `contrib` defined above + coeff = contrib; + } if constexpr (F != FldsID::Nppc) { // for nppc calculation ... // ... do not take volume, weights or smoothing into account @@ -452,6 +430,231 @@ 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] = COORD(i1) + HALF; + x_Code[1] = COORD(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); + + // 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; + } + + 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^ν = 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] }; + + // ZAMO fallback for cells with insufficient particles + if (cmp::AlmostZero(N_norm_sq) || N_norm_sq >= ZERO) { + zamo_fallback_2d(); + 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( + 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] = 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 + 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^ν = 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( + 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] = COORD(i1) + HALF; + x_Code[1] = COORD(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] = 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), + 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/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 4df88cf4c..7670549b1 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -218,6 +218,75 @@ 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 - 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 * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min) * SQR(math::sin(x[1] * dtheta + x2_min)); + } else { + 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 + 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)) { + // 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 + 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_inv * z(x[0] * dr + x1_min, x[1] * dtheta + x2_min); + } else if constexpr (i == 1 && j == 1) { + // g^11 + 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); + } 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 @@ -595,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 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] + norm }; + 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..e1debf689 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 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] + norm }; + 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..beab85b0a 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), theta) * SQR(math::sin(theta)); + } else { + 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 + 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 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] + norm }; + 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 ------------ */ /** 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; +} diff --git a/src/output/fields.cpp b/src/output/fields.cpp index e6b86296f..7e3d05cab 100644 --- a/src/output/fields.cpp +++ b/src/output/fields.cpp @@ -36,27 +36,30 @@ 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: 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 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 = {};