diff --git a/CHANGELOG.md b/CHANGELOG.md index 91254edc3bb3..4ff403842ba9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,12 @@ # Changelog ## Current develop + ### Added (new features/APIs/variables/...) ### Changed (changing behavior/API/variables/...) +- [[PR 1253]](https://github.com/parthenon-hpc-lab/parthenon/pull/1253) Add support for uint64 swarm variables and add default id ### Fixed (not changing behavior/API/variables/...) @@ -17,6 +19,7 @@ ### Incompatibilities (i.e. breaking changes) +- [[PR 1253]](https://github.com/parthenon-hpc-lab/parthenon/pull/1253) Add support for uint64 swarm variables and add default id diff --git a/doc/sphinx/src/interface/metadata.rst b/doc/sphinx/src/interface/metadata.rst index 7b4a76fbba20..4aac247e3bfa 100644 --- a/doc/sphinx/src/interface/metadata.rst +++ b/doc/sphinx/src/interface/metadata.rst @@ -218,6 +218,14 @@ requested. ``FluxRequest::Any`` does not modify search parameters. You will get flux or non-flux variables, and variable associations will be ignored. +``Swarm`` +--------- +By default all partices in a ``Swarm`` contain a persistent, unique identifier +field. If this behavior is not desired (e.g., because anonymous particles are +constantly created and destroyed in a given numerical method), +the ``Metadata::NoPersistentParticleIds`` flag prevent the standard allocation +(and communication of a ``swarm.id`` field). + Application Metadata Flags --------------------------- diff --git a/doc/sphinx/src/particles.rst b/doc/sphinx/src/particles.rst index b0dbbf4944d1..76ab67b284d7 100644 --- a/doc/sphinx/src/particles.rst +++ b/doc/sphinx/src/particles.rst @@ -13,16 +13,20 @@ A ``Swarm`` contains all the particle data for all particles of a given species. It owns a set of ``ParticleVariable``\ s, one for each value of each particle. For example, the spatial positions ``x``, ``y``, and ``z`` of the particles in a swarm are three separate -``ParticleVariable``\ s. ``ParticleVariable``\ s can be either ``Real``- -or ``int``-valued, which is specified by the metadata values -``Metadata::Real`` and ``Metadata::Integer``. ``ParticleVariable``\ s +``ParticleVariable``\ s. ``ParticleVariable``\ s can be either ``Real``-, +``uint64_t``-, or ``int``-valued, which is specified by the metadata values +``Metadata::Real``, ``MetaData::UInt64``, and ``Metadata::Integer``. ``ParticleVariable``\ s should also contain the ``Metadata::Particle`` flag. By default, ``ParticleVariable``\ s provide one scalar quantity per particle, but up to 2D data per particle is currently supported, by passing ``std::vector{N1, N2}`` as the second argument to the ``ParticleVariable`` ``Metadata``. All ``Swarm``\ s by default contain -``x``, ``y``, and ``z`` ``ParticleVariable``\ s; additional fields can -be added as: +positional ``x``, ``y``, and ``z`` ``ParticleVariable``\ s and a ``uint64`` particle +``id`` that is persistent and unique for a given simulation. +The latter field can be disabled (e.g., because it is not required for a method +that constant creates and destroys anonymous particles) by providing the +``Metadata::NoPersistentParticleIds`` flag; +additional fields can be added as: .. code:: cpp @@ -256,17 +260,14 @@ per-swarm list will be output for that swarm. Some visualization tools, like Visit and Paraview, prefer to have access to an ``id`` field for each particle, however it's not clear that a unique ID is required for each particle in - general. Therefore, swarms do not automatically contain an ID swarm - variable. However, when Parthenon outputs a swarm, it automatically - generates an ID variable even if one is not present or - requested. If a variable named ``id'' is available **and** the user - requests it be output, Parthenon will use it. Otherwise, Parthenon - will generate an ``id`` variable just for output and write it to - file. + general. Therefore, if swarms do not automatically contain a unique ID swarm + variable (because the ``Metadata::NoPersistentParticleIds`` has been passed when + creating the swarm), Parthenon will generate an ``id`` variable just for output + (i.e., it is still not available within the simulation itself) and write it to file. .. warning:: - The automatically generted ``id`` is unique for a snapshot in time, + The automatically generated ``id`` is unique for a snapshot in time, but not guaranteed to be time invariant. Indeed it is likely **not** the same between dumps. @@ -280,11 +281,7 @@ Putting it all together, you might have an output block that looks like this: swarms = swarm1, swarm2 swarm_variables = shared_var swarm1_variables = per_swarm_var - swarm2_variables = id The result would be that both ``swarm1`` and ``swarm2`` output the -variables ``x``, ``y``, ``z``, and ``shared_var``. But only ``swarm1`` -outputs ``per_swarm_var``. Both ``swarm1`` and ``swarm2`` will output -an ``id`` field. But the ``id`` field for ``swarm1`` will be -automatically generated, but the ``id`` field for ``swarm2`` will use -the user-initialized value if such a quantity is available. +variables ``id``, ``x``, ``y``, ``z``, and ``shared_var``. But only ``swarm1`` +outputs ``per_swarm_var``. diff --git a/example/particle_leapfrog/particle_leapfrog.cpp b/example/particle_leapfrog/particle_leapfrog.cpp index dc30d03f305a..cfe5c7cf67f6 100644 --- a/example/particle_leapfrog/particle_leapfrog.cpp +++ b/example/particle_leapfrog/particle_leapfrog.cpp @@ -65,7 +65,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { std::string swarm_name = "my_particles"; Metadata swarm_metadata({Metadata::Provides, Metadata::None, Metadata::Independent}); pkg->AddSwarm(swarm_name, swarm_metadata); - pkg->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); Metadata vreal_swarmvalue_metadata({Metadata::Real, Metadata::Vector}, std::vector{3}); pkg->AddSwarmValue("v", swarm_name, vreal_swarmvalue_metadata); @@ -177,7 +176,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto new_particles_context = swarm->AddEmptyParticles(num_particles_this_block); - auto &id = swarm->Get("id").Get(); + auto &id = swarm->Get(swarm_position::id::name()).Get(); auto &x = swarm->Get(swarm_position::x::name()).Get(); auto &y = swarm->Get(swarm_position::y::name()).Get(); auto &z = swarm->Get(swarm_position::z::name()).Get(); @@ -224,16 +223,9 @@ TaskStatus TransportParticles(MeshData *md, const StagedIntegrator *integr swarm_name); auto pack_pos = desc_pos.GetPack(md); - // Make a SwarmPack via strings to get ids - // NOTE(@pdmullen): since we are constructing the pack via strings, we must specify - // the datatype associated with ids (i.e., int). We also extract an indexing map. - std::vector vars_id{"id"}; - static auto desc_id = MakeSwarmPackDescriptor(swarm_name, vars_id); - auto pack_id = desc_id.GetPack(md); - auto pack_id_map = desc_id.GetMap(); - parthenon::PackIdx spi_id(pack_id_map["id"]); - // Make a SwarmPack via strings to get v (note that v is a vector!) + // NOTE(@pdmullen): since we are constructing the pack via strings, we must specify + // the datatype associated with ids (i.e., Real). We also extract an indexing map. std::vector vars_v{"v"}; static auto desc_v = MakeSwarmPackDescriptor(swarm_name, vars_v); auto pack_v = desc_v.GetPack(md); @@ -244,7 +236,6 @@ TaskStatus TransportParticles(MeshData *md, const StagedIntegrator *integr DEFAULT_LOOP_PATTERN, "TestSwarmPack", DevExecSpace(), 0, pack_pos.GetMaxFlatIndex(), KOKKOS_LAMBDA(const int idx) { auto [b, n] = pack_pos.GetBlockParticleIndices(idx); - const int iid = pack_id.GetLowerBound(b, spi_id); const int iv = pack_v.GetLowerBound(b, spi_v); const auto swarm_d = pack_pos.GetContext(b); if (swarm_d.IsActive(n)) { diff --git a/example/particle_tracers/particle_tracers.cpp b/example/particle_tracers/particle_tracers.cpp index d4c10bbee1e6..31996d1ef2c2 100644 --- a/example/particle_tracers/particle_tracers.cpp +++ b/example/particle_tracers/particle_tracers.cpp @@ -30,8 +30,10 @@ #include "bvals/comms/bvals_in_one.hpp" #include "config.hpp" #include "globals.hpp" +#include "interface/metadata.hpp" #include "interface/update.hpp" #include "kokkos_abstraction.hpp" +#include "pack/swarm_default_names.hpp" #include "prolong_restrict/prolong_restrict.hpp" using namespace parthenon::driver::prelude; @@ -146,10 +148,12 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Add swarm of tracer particles std::string swarm_name = "tracers"; - Metadata swarm_metadata({Metadata::Provides, Metadata::None}); + // `NoPersistentParticleIds` is just passed to test this aspect in the regression tests. + // For typical tracers, persistent ids are pretty important. + Metadata swarm_metadata( + {Metadata::Provides, Metadata::None, Metadata::NoPersistentParticleIds}); pkg->AddSwarm(swarm_name, swarm_metadata); Metadata real_swarmvalue_metadata({Metadata::Real}); - pkg->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); pkg->EstimateTimestepBlock = EstimateTimestepBlock; @@ -372,7 +376,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto &x = swarm->Get(swarm_position::x::name()).Get(); auto &y = swarm->Get(swarm_position::y::name()).Get(); auto &z = swarm->Get(swarm_position::z::name()).Get(); - auto &id = swarm->Get("id").Get(); auto swarm_d = swarm->GetDeviceContext(); pmb->par_for( @@ -390,7 +393,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { y(n) = y_min + rng_gen.drand() * (y_max - y_min); z(n) = z_min + rng_gen.drand() * (z_max - z_min); - id(n) = num_tracers * gid + n; rng_pool.free_state(rng_gen); }); diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf.py b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf.py index 85060de24398..7dee275084eb 100644 --- a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf.py +++ b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf.py @@ -57,6 +57,10 @@ def Get(self, variable): print(f"ERROR: {variable} not found in {self.name}") return None + @property + def id(self): + return self.Get("swarm.id") + @property def x(self): return self.Get("swarm.x") diff --git a/src/interface/metadata.hpp b/src/interface/metadata.hpp index 870a34da28fc..0f72cc94fe27 100644 --- a/src/interface/metadata.hpp +++ b/src/interface/metadata.hpp @@ -80,6 +80,8 @@ PARTHENON_INTERNAL_FOR_FLAG(Boolean) \ /** Integer-valued quantity */ \ PARTHENON_INTERNAL_FOR_FLAG(Integer) \ + /** uint64-t-valued quantity */ \ + PARTHENON_INTERNAL_FOR_FLAG(UInt64) \ /** Real-valued quantity */ \ PARTHENON_INTERNAL_FOR_FLAG(Real) \ /************************************************/ \ @@ -123,6 +125,8 @@ /** Align memory of fields to cell centered memory \ (Field will be missing one layer of ghosts if it is not cell centered) **/ \ PARTHENON_INTERNAL_FOR_FLAG(CellMemAligned) \ + /** Particles in a Swarm will not contain a persistent, unique id field **/ \ + PARTHENON_INTERNAL_FOR_FLAG(NoPersistentParticleIds) \ /************************************************/ \ /** Vars specifying coordinates for visualization purposes **/ \ /** You can specify a single 3D var **/ \ @@ -448,6 +452,8 @@ class Metadata { return Boolean; } else if (IsSet(Integer)) { return Integer; + } else if (IsSet(UInt64)) { + return UInt64; } else if (IsSet(Real)) { return Real; } diff --git a/src/interface/swarm.cpp b/src/interface/swarm.cpp index 40632e069716..eb3deb367e89 100644 --- a/src/interface/swarm.cpp +++ b/src/interface/swarm.cpp @@ -11,6 +11,7 @@ // the public, perform publicly and display publicly, and to permit others to do so. //======================================================================================== #include +#include #include #include #include @@ -18,6 +19,7 @@ #include #include +#include "interface/metadata.hpp" #include "mesh/mesh.hpp" #include "pack/swarm_default_names.hpp" #include "swarm.hpp" @@ -84,6 +86,9 @@ Swarm::Swarm(const std::string &label, const Metadata &metadata, const int nmax_ uid_ = get_uid_(label_); // Add default swarm fields + if (!metadata.IsSet(Metadata::NoPersistentParticleIds)) { + Add(swarm_position::id::name(), Metadata({Metadata::UInt64})); + } Add(swarm_position::x::name(), Metadata({Metadata::Real})); Add(swarm_position::y::name(), Metadata({Metadata::Real})); Add(swarm_position::z::name(), Metadata({Metadata::Real})); @@ -118,7 +123,8 @@ void Swarm::Add(const std::string &label, const Metadata &metadata) { // labels must be unique, even between different types of data // if (intMap_.count(label) > 0 || realMap_.count(label) > 0) { if (std::get()>(maps_).count(label) > 0 || - std::get()>(maps_).count(label) > 0) { + std::get()>(maps_).count(label) > 0 || + std::get()>(maps_).count(label) > 0) { throw std::invalid_argument("swarm variable " + label + " already enrolled during Add()!"); } @@ -128,6 +134,8 @@ void Swarm::Add(const std::string &label, const Metadata &metadata) { if (newm.Type() == Metadata::Integer) { Add_(label, newm); + } else if (newm.Type() == Metadata::UInt64) { + Add_(label, newm); } else if (newm.Type() == Metadata::Real) { Add_(label, newm); } else { @@ -145,6 +153,8 @@ void Swarm::Remove(const std::string &label) { auto &int_map = std::get()>(maps_); auto &int_vector = std::get()>(vectors_); + auto &uint64_map = std::get()>(maps_); + auto &uint64_vector = std::get()>(vectors_); auto &real_map = std::get()>(maps_); auto &real_vector = std::get()>(vectors_); @@ -157,7 +167,7 @@ void Swarm::Remove(const std::string &label) { } idx++; } - if (found == true) { + if (found) { // first delete the variable int_vector[idx].reset(); @@ -167,28 +177,44 @@ void Swarm::Remove(const std::string &label) { // Also remove variable from map int_map.erase(label); + return; } - if (found == false) { - idx = 0; - for (const auto &v : real_vector) { - if (label == v->label()) { - found = true; - break; - } - idx++; + // search next variable type (real) + idx = 0; + for (const auto &v : real_vector) { + if (label == v->label()) { + found = true; + break; } + idx++; } - if (found == true) { + if (found) { real_vector[idx].reset(); if (real_vector.size() > 1) real_vector[idx] = std::move(real_vector.back()); real_vector.pop_back(); real_map.erase(label); + return; } - if (found == false) { - throw std::invalid_argument("swarm variable not found in Remove()"); + // search next variable type (uint64_t) + idx = 0; + for (const auto &v : uint64_vector) { + if (label == v->label()) { + found = true; + break; + } + idx++; } + if (found) { + uint64_vector[idx].reset(); + if (uint64_vector.size() > 1) uint64_vector[idx] = std::move(uint64_vector.back()); + uint64_vector.pop_back(); + uint64_map.erase(label); + return; + } + + throw std::invalid_argument("swarm variable not found in Remove()"); } void Swarm::SetPoolMax(const std::int64_t nmax_pool) { @@ -218,6 +244,7 @@ void Swarm::SetPoolMax(const std::int64_t nmax_pool) { pmb->LogMemUsage(n_new * sizeof(int)); auto &int_vector = std::get()>(vectors_); + auto &uint64_vector = std::get()>(vectors_); auto &real_vector = std::get()>(vectors_); for (auto &d : int_vector) { @@ -226,6 +253,12 @@ void Swarm::SetPoolMax(const std::int64_t nmax_pool) { pmb->LogMemUsage(n_new * sizeof(int)); } + for (auto &d : uint64_vector) { + d->data.Resize(d->data.GetDim(6), d->data.GetDim(5), d->data.GetDim(4), + d->data.GetDim(3), d->data.GetDim(2), nmax_pool); + pmb->LogMemUsage(n_new * sizeof(std::uint64_t)); + } + for (auto &d : real_vector) { d->data.Resize(d->data.GetDim(6), d->data.GetDim(5), d->data.GetDim(4), d->data.GetDim(3), d->data.GetDim(2), nmax_pool); @@ -399,17 +432,23 @@ void Swarm::Defrag() { // Get all dynamical variables in swarm auto &int_vector = std::get()>(vectors_); + auto &uint64_vector = std::get()>(vectors_); auto &real_vector = std::get()>(vectors_); PackIndexMap real_imap; PackIndexMap int_imap; + PackIndexMap uint64_imap; auto vreal = PackAllVariables_(real_imap); auto vint = PackAllVariables_(int_imap); + auto vuint64 = PackAllVariables_(uint64_imap); int real_vars_size = real_vector.size(); int int_vars_size = int_vector.size(); + int uint64_vars_size = uint64_vector.size(); auto real_map = real_imap.Map(); auto int_map = int_imap.Map(); + auto uint64_map = uint64_imap.Map(); const int realPackDim = vreal.GetDim(2); const int intPackDim = vint.GetDim(2); + const int uint64PackDim = vuint64.GetDim(2); // Loop over only the active number of particles, and if mask is empty, copy in particle // using address from prefix sum @@ -423,6 +462,9 @@ void Swarm::Defrag() { for (int vidx = 0; vidx < intPackDim; vidx++) { vint(vidx, n) = vint(vidx, nread); } + for (int vidx = 0; vidx < uint64PackDim; vidx++) { + vuint64(vidx, n) = vuint64(vidx, nread); + } mask(n) = true; } }); diff --git a/src/interface/swarm.hpp b/src/interface/swarm.hpp index 49cd906c2a15..3a8ab81072ec 100644 --- a/src/interface/swarm.hpp +++ b/src/interface/swarm.hpp @@ -79,11 +79,14 @@ class Swarm { private: static const int IntVec = 0; static const int RealVec = 1; + static const int UInt64Vec = 2; template static constexpr int getType() { if (std::is_same::value) { return IntVec; + } else if (std::is_same::value) { + return UInt64Vec; } return RealVec; } @@ -214,6 +217,9 @@ class Swarm { for (auto &v : std::get<1>(vectors_)) { size += v->NumComponents(); } + for (auto &v : std::get<2>(vectors_)) { + size += v->NumComponents(); + } return size; } @@ -268,9 +274,11 @@ class Swarm { Metadata m_; int nmax_pool_; std::string info_; - std::tuple, ParticleVariableVector> vectors_; + std::tuple, ParticleVariableVector, + ParticleVariableVector> + vectors_; - std::tuple, MapToParticle> maps_; + std::tuple, MapToParticle, MapToParticle> maps_; ParArray1D mask_; ParArray1D marked_for_removal_; diff --git a/src/interface/swarm_comms.cpp b/src/interface/swarm_comms.cpp index 772eb48161d3..38b7a0c8cbab 100644 --- a/src/interface/swarm_comms.cpp +++ b/src/interface/swarm_comms.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -166,13 +167,17 @@ void Swarm::LoadBuffers_() { pmb->exec_space.fence(); auto &int_vector = std::get()>(vectors_); + auto &uint64_vector = std::get()>(vectors_); auto &real_vector = std::get()>(vectors_); PackIndexMap real_imap; PackIndexMap int_imap; + PackIndexMap uint64_imap; auto vreal = PackAllVariables_(real_imap); auto vint = PackAllVariables_(int_imap); + auto vuint64 = PackAllVariables_(uint64_imap); const int realPackDim = vreal.GetDim(2); const int intPackDim = vint.GetDim(2); + const int uint64PackDim = vuint64.GetDim(2); auto &x = Get(swarm_position::x::name()).Get(); auto &y = Get(swarm_position::y::name()).Get(); @@ -259,10 +264,19 @@ void Swarm::LoadBuffers_() { bdvar.send[bufid](buffer_index) = vreal(i, p_index); buffer_index++; } + // Making sure we catch/update this, when allowing Real = float again + static_assert(sizeof(Real) == 2 * sizeof(int)); for (int i = 0; i < intPackDim; i++) { bdvar.send[bufid](buffer_index) = static_cast(vint(i, p_index)); buffer_index++; } + // Should eventually be a bit_cast once we're on C++20 + static_assert(sizeof(Real) == sizeof(std::uint64_t)); + for (int i = 0; i < uint64PackDim; i++) { + std::memcpy(&bdvar.send[bufid](buffer_index), &vuint64(i, p_index), + sizeof(std::uint64_t)); + buffer_index++; + } } } }); @@ -319,13 +333,17 @@ void Swarm::UnloadBuffers_() { auto neighbor_buffer_index = neighbor_buffer_index_; auto &int_vector = std::get()>(vectors_); + auto &uint64_vector = std::get()>(vectors_); auto &real_vector = std::get()>(vectors_); PackIndexMap real_imap; PackIndexMap int_imap; + PackIndexMap uint64_imap; auto vreal = PackAllVariables_(real_imap); auto vint = PackAllVariables_(int_imap); - int realPackDim = vreal.GetDim(2); - int intPackDim = vint.GetDim(2); + auto vuint64 = PackAllVariables_(uint64_imap); + const int realPackDim = vreal.GetDim(2); + const int intPackDim = vint.GetDim(2); + const int uint64PackDim = vuint64.GetDim(2); const int particle_size = GetParticleDataSize(); auto swarm_d = GetDeviceContext(); @@ -340,10 +358,6 @@ void Swarm::UnloadBuffers_() { } neighbor_received_particles.DeepCopy(neighbor_received_particles_h); - auto &x = Get(swarm_position::x::name()).Get(); - auto &y = Get(swarm_position::y::name()).Get(); - auto &z = Get(swarm_position::z::name()).Get(); - pmb->par_for( PARTHENON_AUTO_LABEL, 0, newParticlesContext.GetNewParticlesMaxIndex(), // n is both new particle index and index over buffer values @@ -363,10 +377,18 @@ void Swarm::UnloadBuffers_() { vreal(i, sid) = bdvar.recv[nbid](bid); bid++; } + // Making sure we catch/update this, when allowing Real = float again + static_assert(sizeof(Real) == 2 * sizeof(int)); for (int i = 0; i < intPackDim; i++) { vint(i, sid) = static_cast(bdvar.recv[nbid](bid)); bid++; } + // Should eventually be a bit_cast once we're on C++20 + static_assert(sizeof(Real) == sizeof(std::uint64_t)); + for (int i = 0; i < uint64PackDim; i++) { + std::memcpy(&vuint64(i, sid), &bdvar.recv[nbid](bid), sizeof(std::uint64_t)); + bid++; + } }); } } diff --git a/src/interface/variable.cpp b/src/interface/variable.cpp index 396edbd73cfe..6050c142249c 100644 --- a/src/interface/variable.cpp +++ b/src/interface/variable.cpp @@ -205,6 +205,7 @@ std::string ParticleVariable::info() const { template class Variable; template class ParticleVariable; template class ParticleVariable; +template class ParticleVariable; template class ParticleVariable; } // namespace parthenon diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp index f5c517cbf5cc..1ad510e2132e 100644 --- a/src/outputs/output_utils.cpp +++ b/src/outputs/output_utils.cpp @@ -187,6 +187,10 @@ AllSwarmInfo::AllSwarmInfo(BlockList_t &block_list, const auto &varname = var->label(); info.Add(varname, var); } + for (const auto &var : swarm->GetVariableVector()) { + const auto &varname = var->label(); + info.Add(varname, var); + } for (const auto &var : swarm->GetVariableVector()) { const auto &varname = var->label(); info.Add(varname, var); @@ -202,6 +206,9 @@ AllSwarmInfo::AllSwarmInfo(BlockList_t &block_list, if (swarm->Contains(varname)) { auto var = swarm->GetP(varname); info.Add(varname, var); + } else if (swarm->Contains(varname)) { + auto var = swarm->GetP(varname); + info.Add(varname, var); } else if (swarm->Contains(varname)) { auto var = swarm->GetP(varname); info.Add(varname, var); diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index b844d7dfbb1a..2f0fbbd4ccab 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -195,7 +195,7 @@ struct SwarmVarInfo { std::array n; int nvar, tensor_rank; bool vector; - std::string swtype; // string for xdmf. "Int" or "Float" + std::string swtype; // string for xdmf. "Int", "UInt" or "Float" SwarmVarInfo() = default; SwarmVarInfo(int n6, int n5, int n4, int n3, int n2, int rank, const std::string &swtype, bool vector) @@ -212,13 +212,14 @@ struct SwarmInfo { SwarmInfo() = default; template using MapToVarVec = std::map>; - std::tuple, MapToVarVec> vars; // SwarmVars on each meshblock - std::map var_info; // size of each swarmvar - std::size_t count_on_rank = 0; // per-meshblock - std::size_t global_offset; // global - std::size_t global_count; // global - std::vector counts; // per-meshblock - std::vector offsets; // global + std::tuple, MapToVarVec, MapToVarVec> + vars; // SwarmVars on each meshblock + std::map var_info; // size of each swarmvar + std::size_t count_on_rank = 0; // per-meshblock + std::size_t global_offset; // global + std::size_t global_count; // global + std::vector counts; // per-meshblock + std::vector offsets; // global // std::vector> masks; // used for reading swarms without defrag std::vector max_indices; // JMM: If we defrag, unneeded? void AddOffsets(const SP_Swarm &swarm); // sets above metadata @@ -233,7 +234,14 @@ struct SwarmInfo { bool vector = m.IsSet(Metadata::Vector); auto shape = m.Shape(); int rank = shape.size(); - std::string t = std::is_same::value ? "Int" : "Float"; + std::string t; + if (std::is_same::value) { + t = "Int"; + } else if (std::is_same::value) { + t = "UInt"; + } else if (std::is_same::value) { + t = "Float"; + } var_info[varname] = SwarmVarInfo(var->GetDim(6), var->GetDim(5), var->GetDim(4), var->GetDim(3), var->GetDim(2), rank, t, vector); } diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index ead3d19dca6b..f72208e7021b 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -256,10 +256,12 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { pin->GetVector(pib->block_name, swname + "_variables"); op.swarms[swname].insert(varnames.begin(), varnames.end()); } - // Always output x, y, and z for swarms so that they work with vis tools. - std::vector coords = {swarm_position::x::name(), - swarm_position::y::name(), - swarm_position::z::name()}; + // Always output id, x, y, and z for swarms so that they work with vis tools. + // Note, it's fine to add the id by default (even though it might not actually + // exist) because only variables that do exists are actually being written. + std::vector coords = { + swarm_position::id::name(), swarm_position::x::name(), + swarm_position::y::name(), swarm_position::z::name()}; op.swarms[swname].insert(coords.begin(), coords.end()); } } @@ -293,8 +295,7 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { #ifdef ENABLE_HDF5 op.write_xdmf = pin->GetOrAddBoolean(op.block_name, "write_xdmf", true); op.write_swarm_xdmf = - (restart) ? false - : pin->GetOrAddBoolean(op.block_name, "write_swarm_xdmf", false); + pin->GetOrAddBoolean(op.block_name, "write_swarm_xdmf", false); pnew_type = new PHDF5Output(op, restart); #else msg << "### FATAL ERROR in Outputs constructor" << std::endl diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 42c28d934702..c4b74a15ad71 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -507,6 +507,14 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm HDF5WriteND(g_var, vname, host_data.data(), vinfo.tensor_rank + 1, local_offset, local_count, global_count, pl_xfer, H5P_DEFAULT); } + auto &uint64_vars = std::get>(swinfo.vars); + for (auto &[vname, swmvarvec] : uint64_vars) { + const auto &vinfo = swinfo.var_info.at(vname); + auto host_data = swinfo.FillHostBuffer(vname, swmvarvec); + SetCounts(swinfo, vinfo); + HDF5WriteND(g_var, vname, host_data.data(), vinfo.tensor_rank + 1, local_offset, + local_count, global_count, pl_xfer, H5P_DEFAULT); + } std::vector pos_tmp; // tmp vector to (potentially) hold particle positions auto &rvars = std::get>(swinfo.vars); for (auto &[vname, swmvarvec] : rvars) { @@ -541,9 +549,9 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm local_count, global_count, pl_xfer, H5P_DEFAULT); } - // If swarm does not contain an "id" object, generate a sequential - // one for vis. - if (swinfo.var_info.count("id") == 0) { + // If swarm does not contain the default id object, generate a sequential + // one for vis called "id" (to differentiate between the default one) + if (swinfo.var_info.count(swarm_position::id::name()) == 0) { std::vector ids(swinfo.global_count); std::iota(std::begin(ids), std::end(ids), swinfo.global_offset); local_offset[0] = swinfo.global_offset; diff --git a/src/outputs/parthenon_xdmf.cpp b/src/outputs/parthenon_xdmf.cpp index c9e534621b1c..fc3a09ae5f5a 100644 --- a/src/outputs/parthenon_xdmf.cpp +++ b/src/outputs/parthenon_xdmf.cpp @@ -249,7 +249,9 @@ void genXDMF(std::string hdfFile, Mesh *pm, SimTime *tm, IndexDomain domain, int ParticleVariableRef(pxdmf, varname, varinfo, swmname, hdfFile, swminfo.global_count); } - if (swminfo.var_info.count("id") == 0) { + // Write info for default (auto-generated) "id"s in case there's not native id field + // for the given swarm. + if (swminfo.var_info.count(swarm_position::id::name()) == 0) { auto swid = SwarmVarInfo(1, 1, 1, 1, 1, 0, "Int", false); ParticleVariableRef(pxdmf, "id", swid, swmname, hdfFile, swminfo.global_count); } @@ -391,7 +393,8 @@ static std::string ParticlePositionRef(const std::string &prefix, const std::string &hdffile, const std::string &datatype, const std::string &extradims, int particle_count) { - std::string precision_string = (datatype == "Float") ? " Precision=\"8\"" : ""; + std::string precision_string = + (datatype == "Float" || datatype == "UInt") ? " Precision=\"8\"" : ""; auto part = StringPrintf("%s\n" @@ -408,15 +411,20 @@ static std::string ParticleDatasetRef(const std::string &prefix, const std::string &hdffile, const std::string &datatype, const std::string &extradims, int particle_count) { - std::string precision_string = (datatype == "Float") ? " Precision=\"8\"" : ""; - auto part = - StringPrintf("%s\n" - "%s %s:/%s/SwarmVars/%s\n" - "%s\n", - prefix.c_str(), extradims.c_str(), particle_count, varname.c_str(), - datatype.c_str(), precision_string.c_str(), prefix.c_str(), - hdffile.c_str(), swmname.c_str(), varname.c_str(), prefix.c_str()); + std::string precision_string = + (datatype == "Float" || datatype == "UInt") ? " Precision=\"8\"" : ""; + // Some vis tools (like Visit) expect an "id" field for particles in the xdmf file. + // If we use the proper internal id field (rather than the autogenerated one), we'll + // override the internal name with the expected name. + auto known_varname = varname == swarm_position::id::name() ? "id" : varname; + auto part = StringPrintf("%s\n" + "%s %s:/%s/SwarmVars/%s\n" + "%s\n", + prefix.c_str(), extradims.c_str(), particle_count, + known_varname.c_str(), datatype.c_str(), + precision_string.c_str(), prefix.c_str(), hdffile.c_str(), + swmname.c_str(), varname.c_str(), prefix.c_str()); return part; } static void ParticleVariableRef(std::ofstream &xdmf, const std::string &varname, diff --git a/src/outputs/restart.hpp b/src/outputs/restart.hpp index 138949dd1667..4f41652fc2eb 100644 --- a/src/outputs/restart.hpp +++ b/src/outputs/restart.hpp @@ -114,6 +114,9 @@ class RestartReader { virtual void ReadSwarmVar(const std::string &swarmname, const std::string &varname, const std::size_t count, const std::size_t offset, const Metadata &m, std::vector &dataVec) = 0; + virtual void ReadSwarmVar(const std::string &swarmname, const std::string &varname, + const std::size_t count, const std::size_t offset, + const Metadata &m, std::vector &dataVec) = 0; virtual void ReadSwarmVar(const std::string &swarmname, const std::string &varname, const std::size_t count, const std::size_t offset, const Metadata &m, std::vector &dataVec) = 0; diff --git a/src/outputs/restart_hdf5.hpp b/src/outputs/restart_hdf5.hpp index 4c346fdac78d..482e17a361a6 100644 --- a/src/outputs/restart_hdf5.hpp +++ b/src/outputs/restart_hdf5.hpp @@ -208,6 +208,11 @@ class RestartReaderHDF5 : public RestartReader { std::vector &dataVec) override { ReadSwarmVar<>(swarmname, varname, count, offset, m, dataVec); }; + void ReadSwarmVar(const std::string &swarmname, const std::string &varname, + const std::size_t count, const std::size_t offset, const Metadata &m, + std::vector &dataVec) override { + ReadSwarmVar<>(swarmname, varname, count, offset, m, dataVec); + }; void ReadSwarmVar(const std::string &swarmname, const std::string &varname, const std::size_t count, const std::size_t offset, const Metadata &m, std::vector &dataVec) override { diff --git a/src/pack/swarm_default_names.hpp b/src/pack/swarm_default_names.hpp index 654da1714550..f272df330184 100644 --- a/src/pack/swarm_default_names.hpp +++ b/src/pack/swarm_default_names.hpp @@ -13,6 +13,7 @@ #ifndef PACK_SWARM_DEFAULT_NAMES_HPP_ #define PACK_SWARM_DEFAULT_NAMES_HPP_ +#include #include #include @@ -28,6 +29,9 @@ namespace swarm_position { +// TODO(review) technically not a position, but it feels overkill to add a separate +// namespace +SWARM_VARIABLE(std::uint64_t, swarm, id); SWARM_VARIABLE(parthenon::Real, swarm, x); SWARM_VARIABLE(parthenon::Real, swarm, y); SWARM_VARIABLE(parthenon::Real, swarm, z); diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 26b5909e8775..e009e1a8ab96 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -397,6 +397,7 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { block_index++; } ReadSwarmVars_(swarm, rm.block_list, count_on_rank, offsets[0]); + ReadSwarmVars_(swarm, rm.block_list, count_on_rank, offsets[0]); ReadSwarmVars_(swarm, rm.block_list, count_on_rank, offsets[0]); } diff --git a/tst/regression/test_suites/particle_leapfrog/parthinput.particle_leapfrog b/tst/regression/test_suites/particle_leapfrog/parthinput.particle_leapfrog index f91bec8010bd..764c2de00208 100644 --- a/tst/regression/test_suites/particle_leapfrog/parthinput.particle_leapfrog +++ b/tst/regression/test_suites/particle_leapfrog/parthinput.particle_leapfrog @@ -56,7 +56,7 @@ cfl = 0.3 file_type = hdf5 dt = 2.0 swarms = my_particles -swarm_variables = id, v, vv +swarm_variables = v, vv file_type = rst diff --git a/tst/regression/test_suites/particle_leapfrog/particle_leapfrog.py b/tst/regression/test_suites/particle_leapfrog/particle_leapfrog.py index fb030a270f06..f90137e274ee 100644 --- a/tst/regression/test_suites/particle_leapfrog/particle_leapfrog.py +++ b/tst/regression/test_suites/particle_leapfrog/particle_leapfrog.py @@ -59,7 +59,7 @@ def Analyse(self, parameters): data = phdf("particles.out0.final.phdf") swarm = data.GetSwarm("my_particles") - inds = np.argsort(swarm["id"]) + inds = np.argsort(swarm.id) final_data = np.vstack((swarm.x, swarm.y, swarm.z, swarm["v"])) final_data = final_data.transpose()[inds] final_data[np.abs(final_data) < 1e-12] = 0 @@ -84,7 +84,11 @@ def Analyse(self, parameters): [0.0, 0.0, 0.0, -1.0, -1.0, -1.0], ] ) + success = True if ref_data.shape != final_data.shape: print("TEST FAIL: Mismatch between actual and reference data shape.") - return False - return (np.abs(final_data - ref_data) <= 1e-10).all() + success = False + if not (np.abs(final_data - ref_data) <= 1e-10).all(): + print("TEST FAIL: Mismatch between actual and reference data shape.") + success = False + return success diff --git a/tst/regression/test_suites/particle_leapfrog_outflow/particle_leapfrog_outflow.py b/tst/regression/test_suites/particle_leapfrog_outflow/particle_leapfrog_outflow.py index 6b855e2a9fbc..8a490a4b1cc5 100644 --- a/tst/regression/test_suites/particle_leapfrog_outflow/particle_leapfrog_outflow.py +++ b/tst/regression/test_suites/particle_leapfrog_outflow/particle_leapfrog_outflow.py @@ -41,7 +41,7 @@ def Analyse(self, parameters): data = phdf("particles.out0.final.phdf") swarm = data.GetSwarm("my_particles") - inds = np.argsort(swarm["id"]) + inds = np.argsort(swarm.id) final_data = np.vstack((swarm.x, swarm.y, swarm.z, swarm["v"])) final_data = final_data.transpose()[inds] diff --git a/tst/regression/test_suites/particle_tracers/particle_tracers.py b/tst/regression/test_suites/particle_tracers/particle_tracers.py index 74cbb123e877..127fbc97a776 100644 --- a/tst/regression/test_suites/particle_tracers/particle_tracers.py +++ b/tst/regression/test_suites/particle_tracers/particle_tracers.py @@ -41,9 +41,11 @@ def Analyse(self, parameters): ) from phdf import phdf + success = True + data = phdf("particle_tracers.out0.final.phdf") swarm = data.GetSwarm("tracers") - inds = np.argsort(swarm["id"]) + inds = np.argsort(swarm.Get("id")) final_data = np.vstack((swarm.x, swarm.y, swarm.z)) final_data = final_data.transpose()[inds] final_data[np.abs(final_data) < 1e-12] = 0 @@ -65,5 +67,16 @@ def Analyse(self, parameters): ) if ref_data.shape != final_data.shape: print("TEST FAIL: Mismatch between actual and reference data shape.") - return False - return (np.abs(final_data - ref_data) <= 1e-8).all() + success = False + if not (np.abs(final_data - ref_data) <= 1e-8).all(): + print( + "Disagreement between reference and actual data:", ref_data, final_data + ) + print( + f"Important: If the numbers match but not the rows, then this is " + f"likely related to a mismatch in the autogereated ids, which are implicitly tested here. " + f"In this case the test here should be updated by ordering particles positions first " + f"and then compare." + ) + success = False + return success