Skip to content

Dev/eckart velocity #187

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Open
alisagk wants to merge 21 commits into
base: 1.4.0rc
Choose a base branch
from
Open

Dev/eckart velocity #187

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
a8b6210
started 4-metric in ks
alisagk Feb 24, 2026
e83ce29
added new functions to support:
alisagk Feb 25, 2026
204d37a
bugs in 4-metric implemetations
alisagk Feb 25, 2026
d9de68e
bug in qks g_ implementation
alisagk Feb 25, 2026
ae8bd68
option for u_0 calculation for photons and massive particles
alisagk Mar 3, 2026
1f81a55
removed dublicate computation of FldsID::V
alisagk Mar 3, 2026
43fbe34
fixes bugs in tensor calculation:
alisagk Mar 3, 2026
88e0f48
removed unnecessary commets
alisagk Mar 3, 2026
e76c03d
stress energy tensor now uses the correct contravariant velocity
alisagk Mar 3, 2026
20a3a05
fixed m^2 and simplified calculation
alisagk Mar 3, 2026
93c3722
added correct u^0 in GR; cut energy definition out -- now uses only u^0
alisagk Mar 3, 2026
fdb99f0
cleaned some variables
alisagk Mar 3, 2026
fa6bfd3
Eckart frame transformation
alisagk Mar 4, 2026
365b33d
Added safety fallback to rest frame Eckart velocity for empty cells
alisagk Mar 4, 2026
f71fbbf
removed unnec calculation
alisagk Mar 9, 2026
d55f170
fixed missing NGHOSTS
alisagk Mar 9, 2026
6e54a61
fallback to ZAMO when no enough particles
alisagk Mar 9, 2026
8974d06
removed unused variable
alisagk Mar 12, 2026
331d641
simplified Tmunu and bulk velocity calculation
alisagk Mar 12, 2026
342c800
always output 4 components of bulk vel in GR
alisagk Mar 12, 2026
36afb1d
removed useless blocks since "V" always outputs 3 components for SR a…
alisagk Mar 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 58 additions & 23 deletions src/framework/domain/output.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -486,18 +486,6 @@ namespace ntt {
{},
local_domain->fields.bckp,
c);
} else if (fld.id() == FldsID::V) {
if constexpr (S != SimEngine::GRPIC) {
ComputeMoments<S, M, FldsID::V>(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);
}
Expand Down Expand Up @@ -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<S, M, FldsID::V>(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<S, M, FldsID::V>(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);
}
Expand Down Expand Up @@ -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<idx_t>(addresses[i]);
raise::ErrorIf(fld.comp[i].size() != 1,
"Wrong # of components requested for 4-velocity",
HERE);
ComputeMoments<S, M, FldsID::V>(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<M::Dim, M, 6>(
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<M::Dim, M, 6>(
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",
Expand Down
3 changes: 2 additions & 1 deletion src/global/global.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,8 @@ namespace ntt {
enum Dimension : unsigned short {
_1D = 1,
_2D = 2,
_3D = 3
_3D = 3,
_4D = 4
};

enum class CellLayer {
Expand Down
Loading
Loading