Skip to content

Commit 69bac9f

Browse files
haykhhaykh
authored
Merge pull request #139 from entity-toolkit/dev/stats
Some issues with the stats output
2 parents 8ff0030 + ba021d6 commit 69bac9f

5 files changed

Lines changed: 107 additions & 73 deletions

File tree

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ set(PROJECT_NAME entity)
77

88
project(
99
${PROJECT_NAME}
10-
VERSION 1.2.3
10+
VERSION 1.2.4
1111
LANGUAGES CXX C)
1212
add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"")
1313
set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ")

‎src/framework/domain/stats.cpp‎

Lines changed: 17 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,7 @@ namespace ntt {
7070
template <SimEngine::type S, class M, StatsID::type P>
7171
auto ComputeMoments(const SimulationParams& params,
7272
const Mesh<M>& mesh,
73+
const M& global_metric,
7374
const std::vector<Particles<M::Dim, M::CoordType>>& prtl_species,
7475
const std::vector<spidx_t>& species,
7576
const std::vector<unsigned short>& components) -> real_t {
@@ -99,6 +100,7 @@ namespace ntt {
99100
if (P == StatsID::Rho and cmp::AlmostZero_host(prtl_spec.mass())) {
100101
continue;
101102
}
103+
real_t temp_buff = ZERO;
102104
Kokkos::parallel_reduce(
103105
"ComputeMoments",
104106
prtl_spec.rangeActiveParticles(),
@@ -111,9 +113,15 @@ namespace ntt {
111113
prtl_spec.mass(), prtl_spec.charge(),
112114
use_weights, mesh.metric),
113115
// clang-format on
114-
buffer);
116+
temp_buff);
117+
buffer += temp_buff;
118+
}
119+
if (P != StatsID::Npart) {
120+
return buffer / (global_metric.totVolume() *
121+
params.template get<real_t>("particles.ppc0"));
122+
} else {
123+
return buffer;
115124
}
116-
return buffer;
117125
}
118126

119127
template <SimEngine::type S, class M, StatsID::type F>
@@ -193,8 +201,8 @@ namespace ntt {
193201
}
194202
auto local_domain = subdomain_ptr(l_subdomain_indices()[0]);
195203
logger::Checkpoint("Writing stats", HERE);
196-
g_stats_writer.write(current_step);
197-
g_stats_writer.write(current_time);
204+
g_stats_writer.write(current_step, false);
205+
g_stats_writer.write(current_time, false);
198206
for (const auto& stat : g_stats_writer.statsWriters()) {
199207
if (stat.id() == StatsID::Custom) {
200208
if (CustomStat != nullptr) {
@@ -206,27 +214,31 @@ namespace ntt {
206214
} else if (stat.id() == StatsID::N) {
207215
g_stats_writer.write(ComputeMoments<S, M, StatsID::N>(params,
208216
local_domain->mesh,
217+
g_mesh.metric,
209218
local_domain->species,
210219
stat.species,
211220
{}));
212221
} else if (stat.id() == StatsID::Npart) {
213222
g_stats_writer.write(
214223
ComputeMoments<S, M, StatsID::Npart>(params,
215224
local_domain->mesh,
225+
g_mesh.metric,
216226
local_domain->species,
217227
stat.species,
218228
{}));
219229
} else if (stat.id() == StatsID::Rho) {
220230
g_stats_writer.write(
221231
ComputeMoments<S, M, StatsID::Rho>(params,
222232
local_domain->mesh,
233+
g_mesh.metric,
223234
local_domain->species,
224235
stat.species,
225236
{}));
226237
} else if (stat.id() == StatsID::Charge) {
227238
g_stats_writer.write(
228239
ComputeMoments<S, M, StatsID::Charge>(params,
229240
local_domain->mesh,
241+
g_mesh.metric,
230242
local_domain->species,
231243
stat.species,
232244
{}));
@@ -235,6 +247,7 @@ namespace ntt {
235247
g_stats_writer.write(
236248
ComputeMoments<S, M, StatsID::T>(params,
237249
local_domain->mesh,
250+
g_mesh.metric,
238251
local_domain->species,
239252
stat.species,
240253
comp));

‎src/kernels/reduced_stats.hpp‎

Lines changed: 74 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -466,82 +466,97 @@ namespace kernel {
466466
}
467467

468468
Inline void operator()(index_t p, real_t& buff) const {
469-
if (tag(p) == ParticleTag::dead) {
469+
if (tag(p) != ParticleTag::alive) {
470470
return;
471471
}
472+
auto dV = ONE;
473+
472474
if constexpr (P == StatsID::Npart) {
473475
buff += ONE;
474476
return;
475-
} else if constexpr (P == StatsID::N or P == StatsID::Rho or
476-
P == StatsID::Charge) {
477-
buff += use_weights ? weight(p) : contrib;
478-
return;
479477
} else {
480-
// for stress-energy tensor
481-
real_t energy { ZERO };
482-
vec_t<Dim::_3D> u_Phys { ZERO };
483-
if constexpr (S == SimEngine::SRPIC) {
484-
// SR
485-
// stress-energy tensor for SR is computed in the tetrad (hatted) basis
486-
if constexpr (M::CoordType == Coord::Cart) {
487-
u_Phys[0] = ux1(p);
488-
u_Phys[1] = ux2(p);
489-
u_Phys[2] = ux3(p);
478+
coord_t<D> x_Code { ZERO };
479+
if constexpr ((D == Dim::_1D) or (D == Dim::_2D) or (D == Dim::_3D)) {
480+
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
481+
}
482+
if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) {
483+
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
484+
}
485+
if constexpr (D == Dim::_3D) {
486+
x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
487+
}
488+
dV = metric.sqrt_det_h(x_Code);
489+
if constexpr (P == StatsID::N or P == StatsID::Rho or P == StatsID::Charge) {
490+
buff += dV * (use_weights ? weight(p) : contrib);
491+
} else {
492+
// for stress-energy tensor
493+
real_t energy { ZERO };
494+
vec_t<Dim::_3D> u_Phys { ZERO };
495+
if constexpr (S == SimEngine::SRPIC) {
496+
// SR
497+
// stress-energy tensor for SR is computed in the tetrad (hatted) basis
498+
if constexpr (M::CoordType == Coord::Cart) {
499+
u_Phys[0] = ux1(p);
500+
u_Phys[1] = ux2(p);
501+
u_Phys[2] = ux3(p);
502+
} else {
503+
static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D");
504+
coord_t<M::PrtlDim> x_Code { ZERO };
505+
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
506+
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
507+
if constexpr (D == Dim::_3D) {
508+
x_Code[2] = static_cast<real_t>(i3(p)) +
509+
static_cast<real_t>(dx3(p));
510+
} else {
511+
x_Code[2] = phi(p);
512+
}
513+
metric.template transform_xyz<Idx::XYZ, Idx::T>(
514+
x_Code,
515+
{ ux1(p), ux2(p), ux3(p) },
516+
u_Phys);
517+
}
518+
if (mass == ZERO) {
519+
energy = NORM(u_Phys[0], u_Phys[1], u_Phys[2]);
520+
} else {
521+
energy = mass * math::sqrt(
522+
ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]));
523+
}
490524
} else {
491-
static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D");
492-
coord_t<M::PrtlDim> x_Code { ZERO };
525+
// GR
526+
// stress-energy tensor for GR is computed in contravariant basis
527+
static_assert(D != Dim::_1D, "GRPIC 1D");
528+
coord_t<D> x_Code { ZERO };
493529
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
494530
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
495531
if constexpr (D == Dim::_3D) {
496532
x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
533+
}
534+
vec_t<Dim::_3D> u_Cntrv { ZERO };
535+
// compute u_i u^i for energy
536+
metric.template transform<Idx::D, Idx::U>(x_Code,
537+
{ ux1(p), ux2(p), ux3(p) },
538+
u_Cntrv);
539+
energy = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) +
540+
u_Cntrv[2] * ux3(p);
541+
if (mass == ZERO) {
542+
energy = math::sqrt(energy);
497543
} else {
498-
x_Code[2] = phi(p);
544+
energy = mass * math::sqrt(ONE + energy);
499545
}
500-
metric.template transform_xyz<Idx::XYZ, Idx::T>(
501-
x_Code,
502-
{ ux1(p), ux2(p), ux3(p) },
503-
u_Phys);
546+
metric.template transform<Idx::U, Idx::PU>(x_Code, u_Cntrv, u_Phys);
504547
}
505-
if (mass == ZERO) {
506-
energy = NORM(u_Phys[0], u_Phys[1], u_Phys[2]);
507-
} else {
508-
energy = mass *
509-
math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]));
510-
}
511-
} else {
512-
// GR
513-
// stress-energy tensor for GR is computed in contravariant basis
514-
static_assert(D != Dim::_1D, "GRPIC 1D");
515-
coord_t<D> x_Code { ZERO };
516-
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
517-
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
518-
if constexpr (D == Dim::_3D) {
519-
x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
520-
}
521-
vec_t<Dim::_3D> u_Cntrv { ZERO };
522-
// compute u_i u^i for energy
523-
metric.template transform<Idx::D, Idx::U>(x_Code,
524-
{ ux1(p), ux2(p), ux3(p) },
525-
u_Cntrv);
526-
energy = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + u_Cntrv[2] * ux3(p);
527-
if (mass == ZERO) {
528-
energy = math::sqrt(energy);
529-
} else {
530-
energy = mass * math::sqrt(ONE + energy);
531-
}
532-
metric.template transform<Idx::U, Idx::PU>(x_Code, u_Cntrv, u_Phys);
533-
}
534-
// compute the corresponding moment
535-
real_t coeff = ONE;
548+
// compute the corresponding moment
549+
real_t coeff = ONE;
536550
#pragma unroll
537-
for (const auto& c : { c1, c2 }) {
538-
if (c == 0) {
539-
coeff *= energy;
540-
} else {
541-
coeff *= u_Phys[c - 1];
551+
for (const auto& c : { c1, c2 }) {
552+
if (c == 0) {
553+
coeff *= energy;
554+
} else {
555+
coeff *= u_Phys[c - 1];
556+
}
542557
}
558+
buff += dV * coeff / energy;
543559
}
544-
buff += coeff / energy;
545560
}
546561
}
547562
};

‎src/output/stats.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,8 @@ namespace stats {
3131
// determine the stats ID
3232
const auto pos = name.find("_");
3333
auto name_raw = (pos == std::string::npos) ? name : name.substr(0, pos);
34-
if ((name_raw[0] != 'E') and (name_raw[0] != 'B') and (name_raw[0] != 'J')) {
34+
if ((name_raw[0] != 'E') and (name_raw[0] != 'B') and
35+
(name_raw[0] != 'J') and (name_raw[0] != 'N')) {
3536
name_raw = name_raw.substr(0, name_raw.find_first_of("0123ijxyzt"));
3637
}
3738
if (StatsID::contains(fmt::toLower(name_raw).c_str())) {

‎src/output/stats.h‎

Lines changed: 13 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -167,17 +167,22 @@ namespace stats {
167167
auto shouldWrite(timestep_t, simtime_t) -> bool;
168168

169169
template <typename T>
170-
inline void write(const T& value) const {
170+
inline void write(const T& value, bool communicate = true) const {
171171
auto tot_value { static_cast<T>(0) };
172172
#if defined(MPI_ENABLED)
173-
MPI_Reduce(&value,
174-
&tot_value,
175-
1,
176-
mpi::get_type<T>(),
177-
MPI_SUM,
178-
MPI_ROOT_RANK,
179-
MPI_COMM_WORLD);
173+
if (communicate) {
174+
MPI_Reduce(&value,
175+
&tot_value,
176+
1,
177+
mpi::get_type<T>(),
178+
MPI_SUM,
179+
MPI_ROOT_RANK,
180+
MPI_COMM_WORLD);
181+
} else {
182+
tot_value = value;
183+
}
180184
#else
185+
(void)communicate;
181186
tot_value = value;
182187
#endif
183188
CallOnce(

0 commit comments

Comments
 (0)