From a60e7fa2fd06ca9bb7e5c7d72f82b034ba1dfc95 Mon Sep 17 00:00:00 2001 From: haykh Date: Mon, 15 Sep 2025 15:10:08 -0700 Subject: [PATCH 1/4] Npart not captured fixed in stats --- src/output/stats.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/output/stats.cpp b/src/output/stats.cpp index 6aa65067..3cd77bfd 100644 --- a/src/output/stats.cpp +++ b/src/output/stats.cpp @@ -31,7 +31,8 @@ namespace stats { // determine the stats ID const auto pos = name.find("_"); auto name_raw = (pos == std::string::npos) ? name : name.substr(0, pos); - if ((name_raw[0] != 'E') and (name_raw[0] != 'B') and (name_raw[0] != 'J')) { + if ((name_raw[0] != 'E') and (name_raw[0] != 'B') and + (name_raw[0] != 'J') and (name_raw[0] != 'N')) { name_raw = name_raw.substr(0, name_raw.find_first_of("0123ijxyzt")); } if (StatsID::contains(fmt::toLower(name_raw).c_str())) { From 48cca9951377294c3743fabad7f901c722f0e8e0 Mon Sep 17 00:00:00 2001 From: haykh Date: Mon, 15 Sep 2025 17:47:38 -0700 Subject: [PATCH 2/4] proper units for prtl moments --- src/framework/domain/stats.cpp | 17 ++++- src/kernels/reduced_stats.hpp | 133 ++++++++++++++++++--------------- 2 files changed, 89 insertions(+), 61 deletions(-) diff --git a/src/framework/domain/stats.cpp b/src/framework/domain/stats.cpp index 6c5bb0ff..0f86554a 100644 --- a/src/framework/domain/stats.cpp +++ b/src/framework/domain/stats.cpp @@ -70,6 +70,7 @@ namespace ntt { template auto ComputeMoments(const SimulationParams& params, const Mesh& mesh, + const M& global_metric, const std::vector>& prtl_species, const std::vector& species, const std::vector& components) -> real_t { @@ -99,6 +100,7 @@ namespace ntt { if (P == StatsID::Rho and cmp::AlmostZero_host(prtl_spec.mass())) { continue; } + real_t temp_buff = ZERO; Kokkos::parallel_reduce( "ComputeMoments", prtl_spec.rangeActiveParticles(), @@ -111,9 +113,15 @@ namespace ntt { prtl_spec.mass(), prtl_spec.charge(), use_weights, mesh.metric), // clang-format on - buffer); + temp_buff); + buffer += temp_buff; + } + if (P != StatsID::Npart) { + return buffer / (global_metric.totVolume() * + params.template get("particles.ppc0")); + } else { + return buffer; } - return buffer; } template @@ -206,6 +214,7 @@ namespace ntt { } else if (stat.id() == StatsID::N) { g_stats_writer.write(ComputeMoments(params, local_domain->mesh, + g_mesh.metric, local_domain->species, stat.species, {})); @@ -213,6 +222,7 @@ namespace ntt { g_stats_writer.write( ComputeMoments(params, local_domain->mesh, + g_mesh.metric, local_domain->species, stat.species, {})); @@ -220,6 +230,7 @@ namespace ntt { g_stats_writer.write( ComputeMoments(params, local_domain->mesh, + g_mesh.metric, local_domain->species, stat.species, {})); @@ -227,6 +238,7 @@ namespace ntt { g_stats_writer.write( ComputeMoments(params, local_domain->mesh, + g_mesh.metric, local_domain->species, stat.species, {})); @@ -235,6 +247,7 @@ namespace ntt { g_stats_writer.write( ComputeMoments(params, local_domain->mesh, + g_mesh.metric, local_domain->species, stat.species, comp)); diff --git a/src/kernels/reduced_stats.hpp b/src/kernels/reduced_stats.hpp index 68e6aa97..fe7f36c5 100644 --- a/src/kernels/reduced_stats.hpp +++ b/src/kernels/reduced_stats.hpp @@ -466,82 +466,97 @@ namespace kernel { } Inline void operator()(index_t p, real_t& buff) const { - if (tag(p) == ParticleTag::dead) { + if (tag(p) != ParticleTag::alive) { return; } + auto dV = ONE; + if constexpr (P == StatsID::Npart) { buff += ONE; return; - } else if constexpr (P == StatsID::N or P == StatsID::Rho or - P == StatsID::Charge) { - buff += use_weights ? weight(p) : contrib; - return; } else { - // for stress-energy tensor - real_t energy { ZERO }; - 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); + coord_t x_Code { ZERO }; + if constexpr ((D == Dim::_1D) or (D == Dim::_2D) or (D == Dim::_3D)) { + x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); + } + if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) { + 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)); + } + dV = metric.sqrt_det_h(x_Code); + if constexpr (P == StatsID::N or P == StatsID::Rho or P == StatsID::Charge) { + buff += dV * (use_weights ? weight(p) : contrib); + } else { + // for stress-energy tensor + real_t energy { ZERO }; + 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 { - static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D"); - coord_t x_Code { ZERO }; + // 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 { - x_Code[2] = phi(p); + energy = mass * math::sqrt(ONE + energy); } - metric.template transform_xyz( - x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); + metric.template transform(x_Code, u_Cntrv, 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 - real_t coeff = ONE; + // compute the corresponding moment + real_t coeff = ONE; #pragma unroll - for (const auto& c : { c1, c2 }) { - if (c == 0) { - coeff *= energy; - } else { - coeff *= u_Phys[c - 1]; + for (const auto& c : { c1, c2 }) { + if (c == 0) { + coeff *= energy; + } else { + coeff *= u_Phys[c - 1]; + } } + buff += dV * coeff / energy; } - buff += coeff / energy; } } }; From ba1374b08da79721908f5940b3871de34a5770c9 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 16 Sep 2025 04:16:53 +0000 Subject: [PATCH 3/4] fixed time/step in stats --- src/framework/domain/stats.cpp | 4 ++-- src/output/stats.h | 20 ++++++++++++-------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/framework/domain/stats.cpp b/src/framework/domain/stats.cpp index 0f86554a..49beadf7 100644 --- a/src/framework/domain/stats.cpp +++ b/src/framework/domain/stats.cpp @@ -201,8 +201,8 @@ namespace ntt { } auto local_domain = subdomain_ptr(l_subdomain_indices()[0]); logger::Checkpoint("Writing stats", HERE); - g_stats_writer.write(current_step); - g_stats_writer.write(current_time); + g_stats_writer.write(current_step, false); + g_stats_writer.write(current_time, false); for (const auto& stat : g_stats_writer.statsWriters()) { if (stat.id() == StatsID::Custom) { if (CustomStat != nullptr) { diff --git a/src/output/stats.h b/src/output/stats.h index fc5bbf3a..5f245cc4 100644 --- a/src/output/stats.h +++ b/src/output/stats.h @@ -167,16 +167,20 @@ namespace stats { auto shouldWrite(timestep_t, simtime_t) -> bool; template - inline void write(const T& value) const { + inline void write(const T& value, bool communicate = true) const { auto tot_value { static_cast(0) }; #if defined(MPI_ENABLED) - MPI_Reduce(&value, - &tot_value, - 1, - mpi::get_type(), - MPI_SUM, - MPI_ROOT_RANK, - MPI_COMM_WORLD); + if (communicate) { + MPI_Reduce(&value, + &tot_value, + 1, + mpi::get_type(), + MPI_SUM, + MPI_ROOT_RANK, + MPI_COMM_WORLD); + } else { + tot_value = value; + } #else tot_value = value; #endif From ba021d660d6a93f5a123277d59cb3543c484fb68 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 16 Sep 2025 09:44:02 -0700 Subject: [PATCH 4/4] version bump + minor warning --- CMakeLists.txt | 2 +- src/output/stats.h | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 29566913..164333ac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ set(PROJECT_NAME entity) project( ${PROJECT_NAME} - VERSION 1.2.3 + VERSION 1.2.4 LANGUAGES CXX C) add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"") set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ") diff --git a/src/output/stats.h b/src/output/stats.h index 5f245cc4..e9161669 100644 --- a/src/output/stats.h +++ b/src/output/stats.h @@ -182,6 +182,7 @@ namespace stats { tot_value = value; } #else + (void)communicate; tot_value = value; #endif CallOnce(