From 1971c1d8ce8f0b8454f7735a013a2b4332712553 Mon Sep 17 00:00:00 2001 From: hayk Date: Tue, 21 Oct 2025 01:55:55 -0400 Subject: [PATCH 01/15] particle checkpoints refactored + added extra vars --- input.example.toml | 13 +- src/checkpoint/reader.cpp | 30 +- src/checkpoint/reader.h | 11 +- src/checkpoint/tests/checkpoint-nompi.cpp | 50 ++- src/checkpoint/writer.cpp | 247 ++++++----- src/checkpoint/writer.h | 30 +- src/framework/CMakeLists.txt | 1 + src/framework/containers/particles.cpp | 31 +- src/framework/containers/particles.h | 53 ++- src/framework/containers/particles_io.cpp | 482 ++++++++++++++++++++++ src/framework/containers/species.h | 50 ++- src/framework/domain/checkpoint.cpp | 337 +-------------- src/framework/parameters.cpp | 27 +- src/framework/tests/parameters.cpp | 25 +- src/framework/tests/particles.cpp | 53 ++- src/output/CMakeLists.txt | 2 + src/output/tests/writer-nompi.cpp | 29 +- src/output/utils/attr_writer.h | 1 + src/output/utils/interpret_prompt.h | 4 +- src/output/utils/readers.cpp | 99 +++++ src/output/utils/readers.h | 48 +++ src/output/utils/writers.cpp | 94 +++++ src/output/utils/writers.h | 55 +++ 23 files changed, 1216 insertions(+), 556 deletions(-) create mode 100644 src/framework/containers/particles_io.cpp create mode 100644 src/output/utils/readers.cpp create mode 100644 src/output/utils/readers.h create mode 100644 src/output/utils/writers.cpp create mode 100644 src/output/utils/writers.h diff --git a/input.example.toml b/input.example.toml index 028102b08..083544647 100644 --- a/input.example.toml +++ b/input.example.toml @@ -302,10 +302,19 @@ # @default: "Boris" [massive]; "Photon" [massless] # @enum: "Boris", "Vay", "Boris,GCA", "Vay,GCA", "Photon", "None" pusher = "" - # Number of additional (payload) variables for each particle of the given species + # Number of additional real-valued variables (payloads) for each particle of the given species # @type: ushort # @default: 0 - n_payloads = "" + n_payloads_real = "" + # Number of additional integer-valued variables (payloads) for each particle of the given species + # @type: ushort + # @default: 0 + # @note: If tracking is enabled, one or two extra integer payloads are reserved (depending on whether MPI is enabled) + n_payloads_int = "" + # Enable tracking of particles using indices for the given species + # @type: bool + # @default: false + tracking = "" # Radiation reaction to use for the species # @type: string # @default: "None" diff --git a/src/checkpoint/reader.cpp b/src/checkpoint/reader.cpp index d973b0ddd..b5d6f0c44 100644 --- a/src/checkpoint/reader.cpp +++ b/src/checkpoint/reader.cpp @@ -108,15 +108,20 @@ namespace checkpoint { } } + template void ReadParticlePayloads(adios2::IO& io, adios2::Engine& reader, + const std::string& suffix, spidx_t s, - array_t& array, + array_t& array, std::size_t nplds, npart_t count, npart_t offset) { - logger::Checkpoint(fmt::format("Reading quantity: s%d_plds", s + 1), HERE); - auto var = io.InquireVariable(fmt::format("s%d_plds", s + 1)); + logger::Checkpoint( + fmt::format("Reading quantity: s%d_pld_%s", s + 1, suffix.c_str()), + HERE); + auto var = io.InquireVariable( + fmt::format("s%d_pld_%s", s + 1, suffix.c_str())); if (var) { var.SetSelection(adios2::Box({ offset, 0 }, { count, nplds })); const auto slice = range_tuple_t(0, count); @@ -126,7 +131,9 @@ namespace checkpoint { adios2::Mode::Sync); Kokkos::deep_copy(array, array_h); } else { - raise::Error(fmt::format("Variable: s%d_plds not found", s + 1), HERE); + raise::Error( + fmt::format("Variable: s%d_pld_%s not found", s + 1, suffix.c_str()), + HERE); } } @@ -158,4 +165,19 @@ namespace checkpoint { CHECKPOINT_PARTICLE_DATA(short) #undef CHECKPOINT_PARTICLE_DATA +#define CHECKPOINT_PARTICLE_PAYLOADS(T) \ + template void ReadParticlePayloads(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + spidx_t, \ + array_t&, \ + std::size_t, \ + npart_t, \ + npart_t); + CHECKPOINT_PARTICLE_PAYLOADS(int) + CHECKPOINT_PARTICLE_PAYLOADS(float) + CHECKPOINT_PARTICLE_PAYLOADS(double) + CHECKPOINT_PARTICLE_PAYLOADS(npart_t) +#undef CHECKPOINT_PARTICLE_DATA + } // namespace checkpoint diff --git a/src/checkpoint/reader.h b/src/checkpoint/reader.h index 7939ba82b..78cc73eb7 100644 --- a/src/checkpoint/reader.h +++ b/src/checkpoint/reader.h @@ -30,11 +30,8 @@ namespace checkpoint { const adios2::Box&, ndfield_t&); - auto ReadParticleCount(adios2::IO&, - adios2::Engine&, - spidx_t, - std::size_t, - std::size_t) -> std::pair; + auto ReadParticleCount(adios2::IO&, adios2::Engine&, spidx_t, std::size_t, std::size_t) + -> std::pair; template void ReadParticleData(adios2::IO&, @@ -45,10 +42,12 @@ namespace checkpoint { npart_t, npart_t); + template void ReadParticlePayloads(adios2::IO&, adios2::Engine&, + const std::string&, spidx_t, - array_t&, + array_t&, std::size_t, npart_t, npart_t); diff --git a/src/checkpoint/tests/checkpoint-nompi.cpp b/src/checkpoint/tests/checkpoint-nompi.cpp index 132a3679a..7be41ca3f 100644 --- a/src/checkpoint/tests/checkpoint-nompi.cpp +++ b/src/checkpoint/tests/checkpoint-nompi.cpp @@ -50,6 +50,11 @@ auto main(int argc, char* argv[]) -> int { array_t i2 { "i_2", npart2 }; array_t u2 { "u_2", npart2 }; + array_t pldr_2 { "pldr_2", npart2, 2 }; + + array_t pldi_1 { "pldi_1", npart1, 1 }; + array_t pldi_2 { "pldi_2", npart2, 2 }; + { // fill data Kokkos::parallel_for( @@ -77,15 +82,20 @@ auto main(int argc, char* argv[]) -> int { "fillPrtl1", npart1, Lambda(index_t p) { - u1(p) = static_cast(p); - i1(p) = static_cast(p); + u1(p) = static_cast(p); + i1(p) = static_cast(p); + pldi_1(p, 0) = static_cast(p * 10); }); Kokkos::parallel_for( "fillPrtl2", npart2, Lambda(index_t p) { - u2(p) = -static_cast(p); - i2(p) = -static_cast(p); + u2(p) = -static_cast(p); + i2(p) = -static_cast(p); + pldr_2(p, 0) = static_cast(p); + pldr_2(p, 1) = static_cast(p * 2); + pldi_2(p, 0) = static_cast(p * 3); + pldi_2(p, 1) = static_cast(p * 4); }); } @@ -101,7 +111,7 @@ auto main(int argc, char* argv[]) -> int { { nx1_gh, nx2_gh, nx3_gh }, { 0, 0, 0 }, { nx1_gh, nx2_gh, nx3_gh }); - writer.defineParticleVariables(Coord::Sph, Dim::_3D, 2, { 0, 2 }); + writer.defineParticleVariables(Coord::Sph, Dim::_3D, 2, { 0, 2 }, { 1, 2 }); writer.beginSaving(0, 0.0); @@ -116,6 +126,11 @@ auto main(int argc, char* argv[]) -> int { writer.saveParticleQuantity("s2_i1", npart2, 0, npart2, i2); writer.saveParticleQuantity("s2_ux1", npart2, 0, npart2, u2); + writer.saveParticlePayloads("s2_pld_r", 2, npart2, 0, npart2, pldr_2); + + writer.saveParticlePayloads("s1_pld_i", 1, npart1, 0, npart1, pldi_1); + writer.saveParticlePayloads("s2_pld_i", 2, npart2, 0, npart2, pldi_2); + writer.endSaving(); } @@ -129,6 +144,11 @@ auto main(int argc, char* argv[]) -> int { array_t i2_read { "i_2", npart2 }; array_t u2_read { "u_2", npart2 }; + array_t pldr_2_read { "pldr_2", npart2, 2 }; + + array_t pldi_1_read { "pldi_1", npart1, 1 }; + array_t pldi_2_read { "pldi_2", npart2, 2 }; + adios2::IO io = adios.DeclareIO("checkpointRead"); adios2::Engine reader = io.Open(checkpoint_path / "step-00000000.bp", adios2::Mode::Read); @@ -147,6 +167,11 @@ auto main(int argc, char* argv[]) -> int { ReadParticleData(io, reader, "i1", 0, i1_read, nprtl1, noff1); ReadParticleData(io, reader, "i1", 1, i2_read, nprtl2, noff2); + ReadParticlePayloads(io, reader, "r", 1, pldr_2_read, 2, nprtl2, noff2); + + ReadParticlePayloads(io, reader, "i", 0, pldi_1_read, 1, nprtl1, noff1); + ReadParticlePayloads(io, reader, "i", 1, pldi_2_read, 2, nprtl2, noff2); + reader.EndStep(); reader.Close(); @@ -182,6 +207,9 @@ auto main(int argc, char* argv[]) -> int { if (i1(p) != i1_read(p)) { raise::KernelError(HERE, "i1 read failed"); } + if (pldi_1(p, 0) != pldi_1_read(p, 0)) { + raise::KernelError(HERE, "pldi_1 read failed"); + } }); Kokkos::parallel_for( "checkPrtl2", @@ -193,6 +221,18 @@ auto main(int argc, char* argv[]) -> int { if (i2(p) != i2_read(p)) { raise::KernelError(HERE, "i2 read failed"); } + if (not cmp::AlmostEqual(pldr_2(p, 0), pldr_2_read(p, 0))) { + raise::KernelError(HERE, "pldr_2(0) read failed"); + } + if (not cmp::AlmostEqual(pldr_2(p, 1), pldr_2_read(p, 1))) { + raise::KernelError(HERE, "pldr_2(1) read failed"); + } + if (pldi_2(p, 0) != pldi_2_read(p, 0)) { + raise::KernelError(HERE, "pldi_2(0) read failed"); + } + if (pldi_2(p, 1) != pldi_2_read(p, 1)) { + raise::KernelError(HERE, "pldi_2(1) read failed"); + } }); } diff --git a/src/checkpoint/writer.cpp b/src/checkpoint/writer.cpp index b766ddfbd..ebb8663ae 100644 --- a/src/checkpoint/writer.cpp +++ b/src/checkpoint/writer.cpp @@ -75,68 +75,84 @@ namespace checkpoint { } } - void Writer::defineParticleVariables(const ntt::Coord& C, - Dimension dim, - std::size_t nspec, - const std::vector& nplds) { - raise::ErrorIf(nplds.size() != nspec, - "Number of payloads does not match the number of species", - HERE); - for (auto s { 0u }; s < nspec; ++s) { - m_io.DefineVariable(fmt::format("s%d_npart", s + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - - for (auto d { 0u }; d < dim; ++d) { - m_io.DefineVariable(fmt::format("s%d_i%d", s + 1, d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - m_io.DefineVariable(fmt::format("s%d_dx%d", s + 1, d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - m_io.DefineVariable(fmt::format("s%d_i%d_prev", s + 1, d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - m_io.DefineVariable(fmt::format("s%d_dx%d_prev", s + 1, d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - } - - if (dim == Dim::_2D and C != ntt::Coord::Cart) { - m_io.DefineVariable(fmt::format("s%d_phi", s + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - } - - for (auto d { 0u }; d < 3; ++d) { - m_io.DefineVariable(fmt::format("s%d_ux%d", s + 1, d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - } - - m_io.DefineVariable(fmt::format("s%d_tag", s + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - m_io.DefineVariable(fmt::format("s%d_weight", s + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - if (nplds[s] > 0) { - m_io.DefineVariable(fmt::format("s%d_plds", s + 1), - { adios2::UnknownDim, nplds[s] }, - { adios2::UnknownDim, 0 }, - { adios2::UnknownDim, nplds[s] }); - } - } - } + // void Writer::defineParticleVariables(const ntt::Coord& C, + // Dimension dim, + // std::size_t nspec, + // const std::vector& npld_r, + // const std::vector& npld_i) { + // raise::ErrorIf( + // npld_r.size() != nspec, + // "Number of real payloads does not match the number of species", + // HERE); + // raise::ErrorIf( + // npld_i.size() != nspec, + // "Number of int payloads does not match the number of species", + // HERE); + // for (auto s { 0u }; s < nspec; ++s) { + // m_io.DefineVariable(fmt::format("s%d_npart", s + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // m_io.DefineVariable(fmt::format("s%d_counter", s + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // + // for (auto d { 0u }; d < dim; ++d) { + // m_io.DefineVariable(fmt::format("s%d_i%d", s + 1, d + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // m_io.DefineVariable(fmt::format("s%d_dx%d", s + 1, d + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // m_io.DefineVariable(fmt::format("s%d_i%d_prev", s + 1, d + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // m_io.DefineVariable(fmt::format("s%d_dx%d_prev", s + 1, d + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // } + // + // if (dim == Dim::_2D and C != ntt::Coord::Cart) { + // m_io.DefineVariable(fmt::format("s%d_phi", s + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // } + // + // for (auto d { 0u }; d < 3; ++d) { + // m_io.DefineVariable(fmt::format("s%d_ux%d", s + 1, d + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // } + // + // m_io.DefineVariable(fmt::format("s%d_tag", s + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // m_io.DefineVariable(fmt::format("s%d_weight", s + 1), + // { adios2::UnknownDim }, + // { adios2::UnknownDim }, + // { adios2::UnknownDim }); + // if (npld_r[s] > 0) { + // m_io.DefineVariable(fmt::format("s%d_pld_r", s + 1), + // { adios2::UnknownDim, npld_r[s] }, + // { adios2::UnknownDim, 0 }, + // { adios2::UnknownDim, npld_r[s] }); + // } + // if (npld_i[s] > 0) { + // m_io.DefineVariable(fmt::format("s%d_pld_i", s + 1), + // { adios2::UnknownDim, npld_i[s] }, + // { adios2::UnknownDim, 0 }, + // { adios2::UnknownDim, npld_i[s] }); + // } + // } + // } auto Writer::shouldSave(timestep_t step, simtime_t time) -> bool { return m_enabled and m_tracker.shouldWrite(step, time); @@ -193,17 +209,6 @@ namespace checkpoint { }); } - template - void Writer::savePerDomainVariable(const std::string& varname, - std::size_t total, - std::size_t offset, - T data) { - auto var = m_io.InquireVariable(varname); - var.SetShape({ total }); - var.SetSelection(adios2::Box({ offset }, { 1 })); - m_writer.Put(var, &data); - } - void Writer::saveAttrs(const ntt::SimulationParams& params, simtime_t time) { CallOnce([&]() { std::ofstream metadata; @@ -228,53 +233,43 @@ namespace checkpoint { adios2::Mode::Sync); } - template - void Writer::saveParticleQuantity(const std::string& quantity, - npart_t glob_total, - npart_t loc_offset, - npart_t loc_size, - const array_t& data) { - const auto slice = range_tuple_t(0, loc_size); - auto var = m_io.InquireVariable(quantity); - - var.SetShape({ glob_total }); - var.SetSelection(adios2::Box({ loc_offset }, { loc_size })); - - auto data_h = Kokkos::create_mirror_view(data); - Kokkos::deep_copy(data_h, data); - auto data_sub = Kokkos::subview(data_h, slice); - m_writer.Put(var, data_sub.data(), adios2::Mode::Sync); - } - - void Writer::saveParticlePayloads(const std::string& quantity, - std::size_t nplds, - npart_t glob_total, - npart_t loc_offset, - npart_t loc_size, - const array_t& data) { - const auto slice = range_tuple_t(0, loc_size); - auto var = m_io.InquireVariable(quantity); - - var.SetShape({ glob_total, nplds }); - var.SetSelection( - adios2::Box({ loc_offset, 0 }, { loc_size, nplds })); - - auto data_h = Kokkos::create_mirror_view(data); - Kokkos::deep_copy(data_h, data); - auto data_sub = Kokkos::subview(data_h, slice, range_tuple_t(0, nplds)); - m_writer.Put(var, data_sub.data(), adios2::Mode::Sync); - } - -#define CHECKPOINT_PERDOMAIN_VARIABLE(T) \ - template void Writer::savePerDomainVariable(const std::string&, \ - std::size_t, \ - std::size_t, \ - T); - CHECKPOINT_PERDOMAIN_VARIABLE(int) - CHECKPOINT_PERDOMAIN_VARIABLE(float) - CHECKPOINT_PERDOMAIN_VARIABLE(double) - CHECKPOINT_PERDOMAIN_VARIABLE(npart_t) -#undef CHECKPOINT_PERDOMAIN_VARIABLE + // template + // void Writer::saveParticleQuantity(const std::string& quantity, + // npart_t glob_total, + // npart_t loc_offset, + // npart_t loc_size, + // const array_t& data) { + // const auto slice = range_tuple_t(0, loc_size); + // auto var = m_io.InquireVariable(quantity); + // + // var.SetShape({ glob_total }); + // var.SetSelection(adios2::Box({ loc_offset }, { loc_size })); + // + // auto data_h = Kokkos::create_mirror_view(data); + // Kokkos::deep_copy(data_h, data); + // auto data_sub = Kokkos::subview(data_h, slice); + // m_writer.Put(var, data_sub.data(), adios2::Mode::Sync); + // } + // + // template + // void Writer::saveParticlePayloads(const std::string& quantity, + // std::size_t nplds, + // npart_t glob_total, + // npart_t loc_offset, + // npart_t loc_size, + // const array_t& data) { + // const auto slice = range_tuple_t(0, loc_size); + // auto var = m_io.InquireVariable(quantity); + // + // var.SetShape({ glob_total, nplds }); + // var.SetSelection( + // adios2::Box({ loc_offset, 0 }, { loc_size, nplds })); + // + // auto data_h = Kokkos::create_mirror_view(data); + // Kokkos::deep_copy(data_h, data); + // auto data_sub = Kokkos::subview(data_h, slice, range_tuple_t(0, nplds)); + // m_writer.Put(var, data_sub.data(), adios2::Mode::Sync); + // } #define CHECKPOINT_FIELD(D, N) \ template void Writer::saveField(const std::string&, \ @@ -287,16 +282,4 @@ namespace checkpoint { CHECKPOINT_FIELD(Dim::_3D, 6) #undef CHECKPOINT_FIELD -#define CHECKPOINT_PARTICLE_QUANTITY(T) \ - template void Writer::saveParticleQuantity(const std::string&, \ - npart_t, \ - npart_t, \ - npart_t, \ - const array_t&); - CHECKPOINT_PARTICLE_QUANTITY(int) - CHECKPOINT_PARTICLE_QUANTITY(float) - CHECKPOINT_PARTICLE_QUANTITY(double) - CHECKPOINT_PARTICLE_QUANTITY(short) -#undef CHECKPOINT_PARTICLE_QUANTITY - } // namespace checkpoint diff --git a/src/checkpoint/writer.h b/src/checkpoint/writer.h index 6f8bc8cb5..2750e2226 100644 --- a/src/checkpoint/writer.h +++ b/src/checkpoint/writer.h @@ -62,35 +62,23 @@ namespace checkpoint { void saveAttrs(const ntt::SimulationParams&, simtime_t); - template - void savePerDomainVariable(const std::string&, std::size_t, std::size_t, T); - template void saveField(const std::string&, const ndfield_t&); - template - void saveParticleQuantity(const std::string&, - npart_t, - npart_t, - npart_t, - const array_t&); - - void saveParticlePayloads(const std::string&, - std::size_t, - npart_t, - npart_t, - npart_t, - const array_t&); - void defineFieldVariables(const ntt::SimEngine&, const std::vector&, const std::vector&, const std::vector&); - void defineParticleVariables(const ntt::Coord&, - Dimension, - std::size_t, - const std::vector&); + [[nodiscard]] + auto io() -> adios2::IO& { + return m_io; + } + + [[nodiscard]] + auto writer() -> adios2::Engine& { + return m_writer; + } [[nodiscard]] auto enabled() const -> bool { diff --git a/src/framework/CMakeLists.txt b/src/framework/CMakeLists.txt index b74d11bec..07328fff6 100644 --- a/src/framework/CMakeLists.txt +++ b/src/framework/CMakeLists.txt @@ -47,6 +47,7 @@ set(SOURCES if(${output}) list(APPEND SOURCES ${SRC_DIR}/domain/output.cpp) list(APPEND SOURCES ${SRC_DIR}/domain/checkpoint.cpp) + list(APPEND SOURCES ${SRC_DIR}/containers/particles_io.cpp) endif() add_library(ntt_framework ${SOURCES}) diff --git a/src/framework/containers/particles.cpp b/src/framework/containers/particles.cpp index d2db9c491..e9e516221 100644 --- a/src/framework/containers/particles.cpp +++ b/src/framework/containers/particles.cpp @@ -22,10 +22,22 @@ namespace ntt { float ch, npart_t maxnpart, const PrtlPusher& pusher, + bool use_tracking, bool use_gca, const Cooling& cooling, - unsigned short npld) - : ParticleSpecies(index, label, m, ch, maxnpart, pusher, use_gca, cooling, npld) { + unsigned short npld_r, + unsigned short npld_i) + : ParticleSpecies(index, + label, + m, + ch, + maxnpart, + pusher, + use_tracking, + use_gca, + cooling, + npld_r, + npld_i) { if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { i1 = array_t { label + "_i1", maxnpart }; @@ -56,8 +68,11 @@ namespace ntt { tag = array_t { label + "_tag", maxnpart }; - if (npld > 0) { - pld = array_t { label + "_pld", maxnpart, npld }; + if (npld_r > 0) { + pld_r = array_t { label + "_pld_r", maxnpart, npld_r }; + } + if (npld_i > 0) { + pld_i = array_t { label + "_pld_i", maxnpart, npld_i }; } if ((D == Dim::_2D) && (C != Coord::Cart)) { @@ -206,8 +221,12 @@ namespace ntt { RemoveDeadInArray(phi, indices_alive); } - if (npld() > 0) { - RemoveDeadInArray(pld, indices_alive); + if (npld_r() > 0) { + RemoveDeadInArray(pld_r, indices_alive); + } + + if (npld_i() > 0) { + RemoveDeadInArray(pld_i, indices_alive); } Kokkos::Experimental::fill( diff --git a/src/framework/containers/particles.h b/src/framework/containers/particles.h index 8ff74be33..c85759035 100644 --- a/src/framework/containers/particles.h +++ b/src/framework/containers/particles.h @@ -23,6 +23,10 @@ #include +#if defined(OUTPUT_ENABLED) + #include +#endif + #include #include @@ -38,6 +42,7 @@ namespace ntt { private: // Number of currently active (used) particles npart_t m_npart { 0 }; + npart_t m_counter { 0 }; bool m_is_sorted { false }; #if !defined(MPI_ENABLED) @@ -60,8 +65,10 @@ namespace ntt { array_t dx1_prev, dx2_prev, dx3_prev; // Array to tag the particles array_t tag; - // Array to store the particle payloads - array_t pld; + // Array to store real-valued payloads + array_t pld_r; + // Array to store integer-valued payloads + array_t pld_i; // phi coordinate (for axisymmetry) array_t phi; @@ -76,9 +83,11 @@ namespace ntt { * @param ch The charge of the species * @param maxnpart The maximum number of allocated particles for the species * @param pusher The pusher assigned for the species + * @param use_tracking Use particle tracking for the species * @param use_gca Use hybrid GCA pusher for the species * @param cooling The cooling mechanism assigned for the species - * @param npld The number of payloads for the species + * @param npld_r The number of real-valued payloads for the species + * @param npld_i The number of integer-valued payloads for the species */ Particles(spidx_t index, const std::string& label, @@ -87,8 +96,10 @@ namespace ntt { npart_t maxnpart, const PrtlPusher& pusher, bool use_gca, + bool use_tracking, const Cooling& cooling, - unsigned short npld = 0); + unsigned short npld_r = 0, + unsigned short npld_i = 0); /** * @brief Constructor for the particle container @@ -102,9 +113,11 @@ namespace ntt { spec.charge(), spec.maxnpart(), spec.pusher(), + spec.use_tracking(), spec.use_gca(), spec.cooling(), - spec.npld()) {} + spec.npld_r(), + spec.npld_i()) {} Particles(const Particles&) = delete; Particles& operator=(const Particles&) = delete; @@ -136,6 +149,17 @@ namespace ntt { return m_npart; } + /** + * @brief Get the particle counter + */ + [[nodiscard]] + auto counter() const -> npart_t { + return m_counter; + } + + /** + * @brief Check if particles are sorted by tag + */ [[nodiscard]] auto is_sorted() const -> bool { return m_is_sorted; @@ -169,8 +193,9 @@ namespace ntt { footprint += sizeof(prtldx_t) * dx2_prev.extent(0); footprint += sizeof(prtldx_t) * dx3_prev.extent(0); footprint += sizeof(short) * tag.extent(0); - footprint += sizeof(real_t) * pld.extent(0) * pld.extent(1); - footprint += sizeof(real_t) * phi.extent(0); + footprint += sizeof(real_t) * pld_r.extent(0) * pld_r.extent(1); + footprint += sizeof(npart_t) * pld_i.extent(0) * pld_i.extent(1); + footprint += sizeof(real_t) * phi.extent(0); return footprint; } @@ -206,6 +231,14 @@ namespace ntt { m_npart = n; } + /** + * @brief Set the particle counter + * @param n The counter value as a npart_t + */ + void set_counter(npart_t n) { + m_counter = n; + } + void set_unsorted() { m_is_sorted = false; } @@ -220,7 +253,11 @@ namespace ntt { */ void SyncHostDevice(); - // void PrintTags(); +#if defined(OUTPUT_ENABLED) + void CheckpointDeclare(adios2::IO&) const; + void CheckpointRead(adios2::IO&, adios2::Engine&, std::size_t, std::size_t); + void CheckpointWrite(adios2::IO&, adios2::Engine&, std::size_t, std::size_t) const; +#endif }; } // namespace ntt diff --git a/src/framework/containers/particles_io.cpp b/src/framework/containers/particles_io.cpp new file mode 100644 index 000000000..870a6e0ae --- /dev/null +++ b/src/framework/containers/particles_io.cpp @@ -0,0 +1,482 @@ +#include "enums.h" +#include "global.h" + +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/log.h" + +#include "framework/containers/particles.h" +#include "output/utils/readers.h" +#include "output/utils/writers.h" + +#include + +#if defined(MPI_ENABLED) + #include +#endif + +namespace ntt { + + template + void Particles::CheckpointDeclare(adios2::IO& io) const { + logger::Checkpoint( + fmt::format("Declaring particle checkpoint for species #%d", index()), + HERE); + io.DefineVariable(fmt::format("s%d_npart", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + io.DefineVariable(fmt::format("s%d_counter", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + for (auto d { 0u }; d < static_cast(D); ++d) { + io.DefineVariable(fmt::format("s%d_i%d", index(), d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + io.DefineVariable(fmt::format("s%d_dx%d", index(), d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + io.DefineVariable(fmt::format("s%d_i%d_prev", index(), d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + io.DefineVariable(fmt::format("s%d_dx%d_prev", index(), d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + + if constexpr (D == Dim::_2D and C != ntt::Coord::Cart) { + io.DefineVariable(fmt::format("s%d_phi", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + + for (auto d { 0u }; d < 3; ++d) { + io.DefineVariable(fmt::format("s%d_ux%d", index(), d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + + io.DefineVariable(fmt::format("s%d_tag", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + io.DefineVariable(fmt::format("s%d_weight", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + if (npld_r() > 0) { + io.DefineVariable(fmt::format("s%d_pld_r", index()), + { adios2::UnknownDim, npld_r() }, + { adios2::UnknownDim, 0 }, + { adios2::UnknownDim, npld_r() }); + } + if (npld_i() > 0) { + io.DefineVariable(fmt::format("s%d_pld_i", index()), + { adios2::UnknownDim, npld_i() }, + { adios2::UnknownDim, 0 }, + { adios2::UnknownDim, npld_i() }); + } + } + + template + void Particles::CheckpointRead(adios2::IO& io, + adios2::Engine& reader, + std::size_t domains_total, + std::size_t domains_offset) { + logger::Checkpoint( + fmt::format("Reading particle checkpoint for species #%d", index()), + HERE); + raise::ErrorIf(npart() > 0, + "Particles already initialized before reading checkpoint", + HERE); + npart_t npart_offset = 0u; + + out::ReadVariable(io, + reader, + fmt::format("s%d_npart", index()), + m_npart, + domains_offset); + + raise::ErrorIf( + npart() > maxnpart(), + fmt::format("npart %d > maxnpart %d after reading checkpoint", + npart(), + maxnpart()), + HERE); + +#if defined(MPI_ENABLED) + { + std::vector glob_nparts(domains_total); + MPI_Allgather(&m_npart, + 1, + mpi::get_type(), + glob_nparts.data(), + 1, + mpi::get_type(), + MPI_COMM_WORLD); + for (auto d { 0u }; d < domains_offset; ++d) { + npart_offset += glob_nparts[d]; + } + } +#endif + out::ReadVariable(io, + reader, + fmt::format("s%d_counter", index()), + m_counter, + domains_offset); + + if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { + out::Read1DArray(io, + reader, + fmt::format("s%d_i1", index()), + i1, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_dx1", index()), + dx1, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_i1_prev", index()), + i1_prev, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_dx1_prev", index()), + dx1_prev, + npart(), + npart_offset); + } + + if constexpr (D == Dim::_2D or D == Dim::_3D) { + out::Read1DArray(io, + reader, + fmt::format("s%d_i2", index()), + i2, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_dx2", index()), + dx2, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_i2_prev", index()), + i2_prev, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_dx2_prev", index()), + dx2_prev, + npart(), + npart_offset); + } + + if constexpr (D == Dim::_3D) { + out::Read1DArray(io, + reader, + fmt::format("s%d_i3", index()), + i3, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_dx3", index()), + dx3, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_i3_prev", index()), + i3_prev, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_dx3_prev", index()), + dx3_prev, + npart(), + npart_offset); + } + + if constexpr (D == Dim::_2D and C != Coord::Cart) { + out::Read1DArray(io, + reader, + fmt::format("s%d_phi", index()), + phi, + npart(), + npart_offset); + } + + out::Read1DArray(io, + reader, + fmt::format("s%d_ux1", index()), + ux1, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_ux2", index()), + ux2, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_ux3", index()), + ux3, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_tag", index()), + tag, + npart(), + npart_offset); + out::Read1DArray(io, + reader, + fmt::format("s%d_weight", index()), + weight, + npart(), + npart_offset); + + if (npld_r() > 0) { + out::Read2DArray(io, + reader, + fmt::format("s%d_pld_r", index()), + pld_r, + npld_r(), + npart(), + npart_offset); + } + + if (npld_i() > 0) { + out::Read2DArray(io, + reader, + fmt::format("s%d_pld_i", index()), + pld_i, + npld_i(), + npart(), + npart_offset); + } + } + + template + void Particles::CheckpointWrite(adios2::IO& io, + adios2::Engine& writer, + std::size_t domains_total, + std::size_t domains_offset) const { + logger::Checkpoint( + fmt::format("Writing particle checkpoint for species #%d", index()), + HERE); + + npart_t npart_offset = 0u; + npart_t npart_total = npart(); + +#if defined(MPI_ENABLED) + { + std::vector glob_nparts(domains_total); + MPI_Allgather(&m_npart, + 1, + mpi::get_type(), + glob_nparts.data(), + 1, + mpi::get_type(), + MPI_COMM_WORLD); + npart_total = 0u; + for (auto r = 0; r < domains_total; ++r) { + if (r < domains_offset) { + npart_offset += glob_nparts[r]; + } + npart_total += glob_nparts[r]; + } + } +#endif + + out::WriteVariable(io, + writer, + fmt::format("s%d_npart", index()), + npart(), + domains_total, + domains_offset); + out::WriteVariable(io, + writer, + fmt::format("s%d_counter", index()), + npart(), + domains_total, + domains_offset); + + if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { + out::Write1DArray(io, + writer, + fmt::format("s%d_i1", index()), + i1, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_dx1", index()), + dx1, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_i1_prev", index()), + i1_prev, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_dx1_prev", index()), + dx1_prev, + npart_total, + npart_offset); + } + + if constexpr (D == Dim::_2D or D == Dim::_3D) { + out::Write1DArray(io, + writer, + fmt::format("s%d_i2", index()), + i2, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_dx2", index()), + dx2, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_i2_prev", index()), + i2_prev, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_dx2_prev", index()), + dx2_prev, + npart_total, + npart_offset); + } + + if constexpr (D == Dim::_3D) { + out::Write1DArray(io, + writer, + fmt::format("s%d_i3", index()), + i3, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_dx3", index()), + dx3, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_i3_prev", index()), + i3_prev, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_dx3_prev", index()), + dx3_prev, + npart_total, + npart_offset); + } + + if constexpr (D == Dim::_2D and C != Coord::Cart) { + out::Write1DArray(io, + writer, + fmt::format("s%d_phi", index()), + phi, + npart_total, + npart_offset); + } + + out::Write1DArray(io, + writer, + fmt::format("s%d_ux1", index()), + ux1, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_ux2", index()), + ux2, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_ux3", index()), + ux3, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_tag", index()), + tag, + npart_total, + npart_offset); + out::Write1DArray(io, + writer, + fmt::format("s%d_weight", index()), + weight, + npart_total, + npart_offset); + if (npld_r() > 0) { + out::Write2DArray(io, + writer, + fmt::format("s%d_pld_r", index()), + pld_r, + npart_total, + npart_offset, + npld_r()); + } + + if (npld_i() > 0) { + out::Write2DArray(io, + writer, + fmt::format("s%d_pld_i", index()), + pld_i, + npart_total, + npart_offset, + npld_i()); + } + } + +#define PARTICLES_CHECKPOINTS(D, C) \ + template void Particles::CheckpointDeclare(adios2::ADIOS&) const; \ + template void Particles::CheckpointRead(adios2::ADIOS&, \ + adios2::Engine&, \ + std::size_t, \ + std::size_t); \ + template void Particles::CheckpointWrite(adios2::IO&, \ + adios2::Engine&, \ + std::size_t, \ + std::size_t) const; \ + PARTICLES_CHECKPOINTS(Dim::_1D, Coord::Cart) \ + PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Cart) \ + PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Cart) \ + PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Sph) \ + PARTICLES_CHECKPOINTS(Dim::_2D, Coord::QSph) \ + PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Sph) \ + PARTICLES_CHECKPOINTS(Dim::_3D, Coord::QSph) +#undef PARTICLES_CHECKPOINTS + +} // namespace ntt diff --git a/src/framework/containers/species.h b/src/framework/containers/species.h index ada0282e2..baf024874 100644 --- a/src/framework/containers/species.h +++ b/src/framework/containers/species.h @@ -33,6 +33,9 @@ namespace ntt { // Pusher assigned for the species const PrtlPusher m_pusher; + // Use particle tracking for the species + const bool m_use_tracking; + // Use byrid gca pusher for the species const bool m_use_gca; @@ -40,7 +43,8 @@ namespace ntt { const Cooling m_cooling; // Number of payloads for the species - const unsigned short m_npld; + const unsigned short m_npld_r; + const unsigned short m_npld_i; public: ParticleSpecies() @@ -50,9 +54,11 @@ namespace ntt { , m_charge { 0.0 } , m_maxnpart { 0 } , m_pusher { PrtlPusher::INVALID } + , m_use_tracking { false } , m_use_gca { false } , m_cooling { Cooling::INVALID } - , m_npld { 0 } {} + , m_npld_r { 0 } + , m_npld_i { 0 } {} /** * @brief Constructor for the particle species container. @@ -63,6 +69,11 @@ namespace ntt { * @param ch The charge of the species. * @param maxnpart The maximum number of allocated particles for the species. * @param pusher The pusher assigned for the species. + * @param use_tracking Use particle tracking for the species. + * @param use_gca Use hybrid GCA pusher for the species. + * @param cooling The cooling mechanism assigned for the species. + * @param npld_r The number of real-valued payloads for the species + * @param npld_i The number of integer-valued payloads for the species */ ParticleSpecies(spidx_t index, const std::string& label, @@ -70,18 +81,35 @@ namespace ntt { float ch, npart_t maxnpart, const PrtlPusher& pusher, + bool use_tracking, bool use_gca, const Cooling& cooling, - unsigned short npld = 0) + unsigned short npld_r = 0, + unsigned short npld_i = 0) : m_index { index } , m_label { std::move(label) } , m_mass { m } , m_charge { ch } , m_maxnpart { maxnpart } , m_pusher { pusher } + , m_use_tracking { use_tracking } , m_use_gca { use_gca } , m_cooling { cooling } - , m_npld { npld } {} + , m_npld_r { npld_r } + , m_npld_i { npld_i } { + if (use_tracking) { +#if !defined(MPI_ENABLED) + raise::ErrorIf(m_npld_i < 1, + "npld_i must be at least 1 when tracking is enabled", + HERE); +#else + raise::ErrorIf( + m_npld_i < 2, + "npld_i must be at least 2 when tracking is enabled with MPI", + HERE); +#endif + } + } ParticleSpecies(const ParticleSpecies&) = default; @@ -120,6 +148,11 @@ namespace ntt { return m_pusher; } + [[nodiscard]] + auto use_tracking() const -> bool { + return m_use_tracking; + } + [[nodiscard]] auto use_gca() const -> bool { return m_use_gca; @@ -131,8 +164,13 @@ namespace ntt { } [[nodiscard]] - auto npld() const -> unsigned short { - return m_npld; + auto npld_r() const -> unsigned short { + return m_npld_r; + } + + [[nodiscard]] + auto npld_i() const -> unsigned short { + return m_npld_i; } }; } // namespace ntt diff --git a/src/framework/domain/checkpoint.cpp b/src/framework/domain/checkpoint.cpp index e0e34f993..8fcefa989 100644 --- a/src/framework/domain/checkpoint.cpp +++ b/src/framework/domain/checkpoint.cpp @@ -42,9 +42,10 @@ namespace ntt { } auto loc_shape_with_ghosts = local_domain->mesh.n_all(); - std::vector nplds; + std::vector npld_r, npld_i; for (auto s { 0u }; s < local_domain->species.size(); ++s) { - nplds.push_back(local_domain->species[s].npld()); + npld_r.push_back(local_domain->species[s].npld_r()); + npld_i.push_back(local_domain->species[s].npld_i()); } const path_t checkpoint_root = params.template get( @@ -62,10 +63,9 @@ namespace ntt { glob_shape_with_ghosts, off_ncells_with_ghosts, loc_shape_with_ghosts); - g_checkpoint_writer.defineParticleVariables(M::CoordType, - M::Dim, - local_domain->species.size(), - nplds); + for (auto& species : local_domain->species) { + species.CheckpointDeclare(g_checkpoint_writer.io()); + } } } @@ -96,165 +96,17 @@ namespace ntt { g_checkpoint_writer.saveField("em0", local_domain->fields.em0); g_checkpoint_writer.saveField("cur0", local_domain->fields.cur0); } - std::size_t dom_offset = 0, dom_tot = 1; + std::size_t dom_tot = 1, dom_offset = 0; #if defined(MPI_ENABLED) - dom_offset = g_mpi_rank; dom_tot = g_mpi_size; + dom_offset = g_mpi_rank; #endif // MPI_ENABLED - for (auto s { 0u }; s < local_domain->species.size(); ++s) { - auto npart = local_domain->species[s].npart(); - npart_t offset = 0; - auto glob_tot = npart; -#if defined(MPI_ENABLED) - auto glob_npart = std::vector(g_ndomains); - MPI_Allgather(&npart, - 1, - mpi::get_type(), - glob_npart.data(), - 1, - mpi::get_type(), - MPI_COMM_WORLD); - glob_tot = 0; - for (auto r = 0; r < g_mpi_size; ++r) { - if (r < g_mpi_rank) { - offset += glob_npart[r]; - } - glob_tot += glob_npart[r]; - } -#endif // MPI_ENABLED - g_checkpoint_writer.savePerDomainVariable( - fmt::format("s%d_npart", s + 1), - dom_tot, - dom_offset, - npart); - if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or - M::Dim == Dim::_3D) { - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_i1", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].i1); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_dx1", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].dx1); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_i1_prev", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].i1_prev); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_dx1_prev", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].dx1_prev); - } - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_i2", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].i2); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_dx2", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].dx2); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_i2_prev", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].i2_prev); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_dx2_prev", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].dx2_prev); - } - if constexpr (M::Dim == Dim::_3D) { - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_i3", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].i3); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_dx3", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].dx3); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_i3_prev", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].i3_prev); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_dx3_prev", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].dx3_prev); - } - if constexpr (M::Dim == Dim::_2D and M::CoordType != Coord::Cart) { - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_phi", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].phi); - } - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_ux1", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].ux1); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_ux2", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].ux2); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_ux3", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].ux3); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_tag", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].tag); - g_checkpoint_writer.saveParticleQuantity( - fmt::format("s%d_weight", s + 1), - glob_tot, - offset, - npart, - local_domain->species[s].weight); - - auto nplds = local_domain->species[s].npld(); - if (nplds > 0) { - g_checkpoint_writer.saveParticlePayloads(fmt::format("s%d_plds", s + 1), - nplds, - glob_tot, - offset, - npart, - local_domain->species[s].pld); - } + for (const auto& species : local_domain->species) { + species.CheckpointWrite(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + dom_tot, + dom_offset); } } g_checkpoint_writer.endSaving(); @@ -284,8 +136,8 @@ namespace ntt { #endif reader.BeginStep(); - for (auto& ldidx : l_subdomain_indices()) { - auto& domain = g_subdomains[ldidx]; + for (const auto local_domain_idx : l_subdomain_indices()) { + auto& domain = g_subdomains[local_domain_idx]; adios2::Box range; for (auto d { 0u }; d < M::Dim; ++d) { range.first.push_back(domain.offset_ncells()[d] + @@ -315,163 +167,10 @@ namespace ntt { range3, domain.fields.cur0); } - for (auto s { 0u }; s < domain.species.size(); ++s) { - const auto [loc_npart, offset_npart] = - checkpoint::ReadParticleCount(io, reader, s, ldidx, ndomains()); - raise::ErrorIf(loc_npart > domain.species[s].maxnpart(), - "loc_npart > domain.species[s].maxnpart()", - HERE); - if (loc_npart == 0) { - continue; - } - if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or - M::Dim == Dim::_3D) { - checkpoint::ReadParticleData(io, - reader, - "i1", - s, - domain.species[s].i1, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "dx1", - s, - domain.species[s].dx1, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "i1_prev", - s, - domain.species[s].i1_prev, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "dx1_prev", - s, - domain.species[s].dx1_prev, - loc_npart, - offset_npart); - } - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - checkpoint::ReadParticleData(io, - reader, - "i2", - s, - domain.species[s].i2, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "dx2", - s, - domain.species[s].dx2, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "i2_prev", - s, - domain.species[s].i2_prev, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "dx2_prev", - s, - domain.species[s].dx2_prev, - loc_npart, - offset_npart); - } - if constexpr (M::Dim == Dim::_3D) { - checkpoint::ReadParticleData(io, - reader, - "i3", - s, - domain.species[s].i3, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "dx3", - s, - domain.species[s].dx3, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "i3_prev", - s, - domain.species[s].i3_prev, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "dx3_prev", - s, - domain.species[s].dx3_prev, - loc_npart, - offset_npart); - } - if constexpr (M::Dim == Dim::_2D and M::CoordType != Coord::Cart) { - checkpoint::ReadParticleData(io, - reader, - "phi", - s, - domain.species[s].phi, - loc_npart, - offset_npart); - } - checkpoint::ReadParticleData(io, - reader, - "ux1", - s, - domain.species[s].ux1, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "ux2", - s, - domain.species[s].ux2, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "ux3", - s, - domain.species[s].ux3, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "tag", - s, - domain.species[s].tag, - loc_npart, - offset_npart); - checkpoint::ReadParticleData(io, - reader, - "weight", - s, - domain.species[s].weight, - loc_npart, - offset_npart); - const auto nplds = domain.species[s].npld(); - if (nplds > 0) { - checkpoint::ReadParticlePayloads(io, - reader, - s, - domain.species[s].pld, - nplds, - loc_npart, - offset_npart); - } - domain.species[s].set_npart(loc_npart); - } // species loop + for (auto& species : domain.species) { + species.CheckpointRead(io, reader, local_domain_idx, ndomains()); + } } // local subdomain loop diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 884551e2c..27906b6e4 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -200,10 +200,21 @@ namespace ntt { const auto maxnpart_real = toml::find(sp, "maxnpart"); const auto maxnpart = static_cast(maxnpart_real); auto pusher = toml::find_or(sp, "pusher", std::string(def_pusher)); - const auto npayloads = toml::find_or(sp, - "n_payloads", - static_cast(0)); - const auto cooling = toml::find_or(sp, "cooling", std::string("None")); + const auto npayloads_real = toml::find_or(sp, + "n_payloads_real", + static_cast(0)); + const auto use_tracking = toml::find_or(sp, "tracking", false); + auto npayloads_int = toml::find_or(sp, + "n_payloads_int", + static_cast(0)); + if (use_tracking) { +#if !defined(MPI_ENABLED) + npayloads_int += 1; +#else + npayloads_int += 2; +#endif + } + const auto cooling = toml::find_or(sp, "cooling", std::string("None")); raise::ErrorIf((fmt::toLower(cooling) != "none") && is_massless, "cooling is only applicable to massive particles", HERE); @@ -241,9 +252,11 @@ namespace ntt { charge, maxnpart, pusher_enum, + use_tracking, use_gca, cooling_enum, - npayloads)); + npayloads_real, + npayloads_int)); idx += 1; } set("particles.species", species); @@ -468,9 +481,11 @@ namespace ntt { particle_species.charge(), maxnpart, particle_species.pusher(), + particle_species.use_tracking(), particle_species.use_gca(), particle_species.cooling(), - particle_species.npld()); + particle_species.npld_r(), + particle_species.npld_i()); idxM1++; } set("particles.species", new_species); diff --git a/src/framework/tests/parameters.cpp b/src/framework/tests/parameters.cpp index 07b2c11b3..c9475d688 100644 --- a/src/framework/tests/parameters.cpp +++ b/src/framework/tests/parameters.cpp @@ -12,7 +12,6 @@ #include #include -#include using namespace toml::literals::toml_literals; const auto mink_1d = u8R"( @@ -55,7 +54,8 @@ const auto mink_1d = u8R"( charge = -1.0 maxnpart = 1e2 pusher = "boris" - n_payloads = 3 + n_payloads_real = 3 + tracking = true [[particles.species]] label = "p+" @@ -139,7 +139,7 @@ const auto sph_2d = u8R"( charge = -1.0 maxnpart = 1e2 pusher = "boris,gca" - n_payloads = 3 + n_payloads_real = 3 cooling = "synchrotron" [[particles.species]] @@ -149,6 +149,7 @@ const auto sph_2d = u8R"( maxnpart = 1e2 pusher = "boris,gca" cooling = "synchrotron" + n_payloads_int = 2 [[particles.species]] label = "ph" @@ -297,7 +298,9 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[0].pusher(), PrtlPusher::BORIS, "species[0].pusher"); - assert_equal(species[0].npld(), 3, "species[0].npld"); + assert_equal(species[0].npld_r(), 3, "species[0].npld_r"); + assert_equal(species[0].npld_i(), 1, "species[0].npld_i"); + assert_equal(species[0].use_tracking(), true, "species[0].tracking"); assert_equal(species[1].label(), "p+", "species[1].label"); assert_equal(species[1].mass(), 1.0f, "species[1].mass"); @@ -306,7 +309,7 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[1].pusher(), PrtlPusher::VAY, "species[1].pusher"); - assert_equal(species[1].npld(), 0, "species[1].npld"); + assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); assert_equal(params_mink_1d.get("setup.myfloat"), (real_t)(1e-2), @@ -417,7 +420,7 @@ auto main(int argc, char* argv[]) -> int { PrtlPusher::BORIS, "species[0].pusher"); assert_equal(species[0].use_gca(), true, "species[0].use_gca"); - assert_equal(species[0].npld(), 3, "species[0].npld"); + assert_equal(species[0].npld_r(), 3, "species[0].npld_r"); assert_equal(species[0].cooling(), Cooling::SYNCHROTRON, "species[0].cooling"); @@ -430,10 +433,12 @@ auto main(int argc, char* argv[]) -> int { PrtlPusher::BORIS, "species[1].pusher"); assert_equal(species[1].use_gca(), true, "species[1].use_gca"); - assert_equal(species[1].npld(), 0, "species[1].npld"); + assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); assert_equal(species[1].cooling(), Cooling::SYNCHROTRON, "species[1].cooling"); + assert_equal(species[1].npld_i(), 2, "species[1].npld_i"); + assert_equal(species[1].use_tracking(), false, "species[1].tracking"); assert_equal(species[2].label(), "ph", "species[2].label"); assert_equal(species[2].mass(), 0.0f, "species[2].mass"); @@ -442,7 +447,7 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[2].pusher(), PrtlPusher::PHOTON, "species[2].pusher"); - assert_equal(species[2].npld(), 0, "species[2].npld"); + assert_equal(species[2].npld_r(), 0, "species[2].npld_r"); } { @@ -551,7 +556,7 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[0].pusher(), PrtlPusher::BORIS, "species[0].pusher"); - assert_equal(species[0].npld(), 0, "species[0].npld"); + assert_equal(species[0].npld_r(), 0, "species[0].npld_r"); assert_equal(species[1].label(), "e+", "species[1].label"); assert_equal(species[1].mass(), 1.0f, "species[1].mass"); @@ -560,7 +565,7 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[1].pusher(), PrtlPusher::BORIS, "species[1].pusher"); - assert_equal(species[1].npld(), 0, "species[1].npld"); + assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); } } catch (std::exception& err) { diff --git a/src/framework/tests/particles.cpp b/src/framework/tests/particles.cpp index 6c4c227b5..ab3aa0e04 100644 --- a/src/framework/tests/particles.cpp +++ b/src/framework/tests/particles.cpp @@ -4,24 +4,33 @@ #include "global.h" #include "utils/error.h" -#include "utils/formatting.h" #include -#include #include -#include template -void testParticles(const int& index, +void testParticles(int index, const std::string& label, - const float& m, - const float& ch, - const std::size_t& maxnpart, + float m, + float ch, + std::size_t maxnpart, const ntt::PrtlPusher& pusher, + bool use_tracking, const ntt::Cooling& cooling, - const unsigned short& npld = 0) { + unsigned short npld_r = 0, + unsigned short npld_i = 0) { using namespace ntt; - auto p = Particles(index, label, m, ch, maxnpart, pusher, false, cooling, npld); + auto p = Particles(index, + label, + m, + ch, + maxnpart, + pusher, + use_tracking, + false, + cooling, + npld_r, + npld_i); raise::ErrorIf(p.index() != index, "Index mismatch", HERE); raise::ErrorIf(p.label() != label, "Label mismatch", HERE); raise::ErrorIf(p.mass() != m, "Mass mismatch", HERE); @@ -46,9 +55,13 @@ void testParticles(const int& index, raise::ErrorIf(p.tag.extent(0) != maxnpart, "tag incorrectly allocated", HERE); raise::ErrorIf(p.weight.extent(0) != maxnpart, "weight incorrectly allocated", HERE); - if (npld > 0) { - raise::ErrorIf(p.pld.extent(0) != maxnpart, "pld incorrectly allocated", HERE); - raise::ErrorIf(p.pld.extent(1) != npld, "pld incorrectly allocated", HERE); + if (npld_r > 0) { + raise::ErrorIf(p.pld_r.extent(0) != maxnpart, "pld_r incorrectly allocated", HERE); + raise::ErrorIf(p.pld_r.extent(1) != npld_r, "pld_r incorrectly allocated", HERE); + } + if (npld_i > 0) { + raise::ErrorIf(p.pld_i.extent(0) != maxnpart, "pld_i incorrectly allocated", HERE); + raise::ErrorIf(p.pld_i.extent(1) != npld_i, "pld_i incorrectly allocated", HERE); } if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) { @@ -103,6 +116,7 @@ auto main(int argc, char** argv) -> int { -1.0, 100, PrtlPusher::BORIS, + false, Cooling::SYNCHROTRON); testParticles(2, "p+", @@ -110,13 +124,17 @@ auto main(int argc, char** argv) -> int { -1.0, 1000, PrtlPusher::VAY, - Cooling::SYNCHROTRON); + true, + Cooling::SYNCHROTRON, + 2, + 1); testParticles(3, "ph", 0.0, 0.0, 100, PrtlPusher::PHOTON, + false, Cooling::NONE, 5); testParticles(4, @@ -125,15 +143,20 @@ auto main(int argc, char** argv) -> int { 1.0, 100, PrtlPusher::BORIS, - Cooling::NONE); + true, + Cooling::NONE, + 2, + 3); testParticles(5, "e+", 1.0, 1.0, 100, PrtlPusher::BORIS, + false, Cooling::NONE, - 1); + 1, + 2); } catch (const std::exception& e) { std::cerr << "Error: " << e.what() << std::endl; Kokkos::finalize(); diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt index 1b132fb60..ef20dae56 100644 --- a/src/output/CMakeLists.txt +++ b/src/output/CMakeLists.txt @@ -30,6 +30,8 @@ set(SOURCES ${SRC_DIR}/stats.cpp ${SRC_DIR}/fields.cpp ${SRC_DIR}/utils/interpret_prompt.cpp) if(${output}) list(APPEND SOURCES ${SRC_DIR}/writer.cpp) + list(APPEND SOURCES ${SRC_DIR}/utils/writers.cpp) + list(APPEND SOURCES ${SRC_DIR}/utils/readers.cpp) endif() add_library(ntt_output ${SOURCES}) diff --git a/src/output/tests/writer-nompi.cpp b/src/output/tests/writer-nompi.cpp index 66d834f43..a185b1777 100644 --- a/src/output/tests/writer-nompi.cpp +++ b/src/output/tests/writer-nompi.cpp @@ -18,8 +18,8 @@ using namespace ntt; void cleanup() { namespace fs = std::filesystem; - fs::path tempfile_path { "test.h5" }; - fs::remove(tempfile_path); + fs::path tempfile_path { "test.bp" }; + fs::remove_all(tempfile_path); } #define CEILDIV(a, b) \ @@ -71,7 +71,7 @@ auto main(int argc, char* argv[]) -> int { { // write auto writer = out::Writer(); - writer.init(&adios, "hdf5", "test", false); + writer.init(&adios, "bpfile", "test", false); writer.defineMeshLayout({ nx1, nx2, nx3 }, { 0, 0, 0 }, { nx1, nx2, nx3 }, @@ -100,19 +100,20 @@ auto main(int argc, char* argv[]) -> int { { // read adios2::IO io = adios.DeclareIO("read-test"); - io.SetEngine("hdf5"); - adios2::Engine reader = io.Open("test.h5", adios2::Mode::Read); - const auto layoutRight = io.InquireAttribute("LayoutRight").Data()[0] == - 1; - - raise::ErrorIf(io.InquireAttribute("NGhosts").Data()[0] != 0, - "NGhosts is not correct", - HERE); - raise::ErrorIf(io.InquireAttribute("Dimension").Data()[0] != 3, - "Dimension is not correct", - HERE); + io.SetEngine("BPFile"); + adios2::Engine reader = io.Open("test.bp", adios2::Mode::Read); for (auto step = 0u; reader.BeginStep() == adios2::StepStatus::OK; ++step) { + const auto layoutRight = io.InquireAttribute("LayoutRight").Data()[0] == + 1; + + raise::ErrorIf(io.InquireAttribute("NGhosts").Data()[0] != 0, + "NGhosts is not correct", + HERE); + raise::ErrorIf(io.InquireAttribute("Dimension").Data()[0] != 3, + "Dimension is not correct", + HERE); + timestep_t step_read; simtime_t time_read; diff --git a/src/output/utils/attr_writer.h b/src/output/utils/attr_writer.h index c8b21e4c2..41b64fb75 100644 --- a/src/output/utils/attr_writer.h +++ b/src/output/utils/attr_writer.h @@ -7,6 +7,7 @@ * @namespaces: * - out:: */ + #ifndef OUTPUT_UTILS_ATTR_WRITER_H #define OUTPUT_UTILS_ATTR_WRITER_H diff --git a/src/output/utils/interpret_prompt.h b/src/output/utils/interpret_prompt.h index 032482cf8..488d81101 100644 --- a/src/output/utils/interpret_prompt.h +++ b/src/output/utils/interpret_prompt.h @@ -26,8 +26,8 @@ namespace out { auto InterpretSpecies(const std::string&) -> std::vector; - auto InterpretComponents( - const std::vector&) -> std::vector>; + auto InterpretComponents(const std::vector&) + -> std::vector>; } // namespace out diff --git a/src/output/utils/readers.cpp b/src/output/utils/readers.cpp new file mode 100644 index 000000000..d2e866af1 --- /dev/null +++ b/src/output/utils/readers.cpp @@ -0,0 +1,99 @@ +#include "output/utils/readers.h" + +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/error.h" +#include "utils/formatting.h" + +#include + +#include + +namespace out { + + template + void ReadVariable(adios2::IO& io, + adios2::Engine& reader, + const std::string& quantity, + T& data, + std::size_t local_offset) { + auto var = io.InquireVariable(quantity); + if (var) { + var.SetSelection(adios2::Box({ local_offset }, { 1 })); + reader.Get(var, &data, adios2::Mode::Sync); + } else { + raise::Error(fmt::format("Variable: %s not found", quantity.c_str()), HERE); + } + } + + template + void Read1DArray(adios2::IO& io, + adios2::Engine& reader, + const std::string& quantity, + array_t& data, + std::size_t local_size, + std::size_t local_offset) { + auto var = io.InquireVariable(quantity); + if (var) { + var.SetSelection(adios2::Box({ local_offset }, { local_size })); + const auto slice = range_tuple_t(0, local_size); + auto data_h = Kokkos::create_mirror_view(data); + reader.Get(var, Kokkos::subview(data_h, slice).data(), adios2::Mode::Sync); + Kokkos::deep_copy(Kokkos::subview(data, slice), + Kokkos::subview(data_h, slice)); + } else { + raise::Error(fmt::format("Variable: %s not found", quantity.c_str()), HERE); + } + } + + template + void Read2DArray(adios2::IO& io, + adios2::Engine& reader, + const std::string& quantity, + array_t& data, + unsigned short dim2_size, + std::size_t local_size, + std::size_t local_offset) { + auto var = io.InquireVariable(quantity); + if (var) { + var.SetSelection(adios2::Box({ local_offset, 0 }, + { local_size, dim2_size })); + const auto slice = range_tuple_t(0, local_size); + auto data_h = Kokkos::create_mirror_view(data); + reader.Get(var, + Kokkos::subview(data_h, slice, range_tuple_t(0, dim2_size)).data(), + adios2::Mode::Sync); + Kokkos::deep_copy(data, data_h); + } else { + raise::Error(fmt::format("Variable: %s not found", quantity.c_str()), HERE); + } + } + +#define ARRAY_READERS(T) \ + template void ReadVariable(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + T&, \ + std::size_t); \ + template void Read1DArray(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + array_t&, \ + std::size_t, \ + std::size_t); \ + template void Read2DArray(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + array_t&, \ + unsigned short, \ + std::size_t, \ + std::size_t); \ + ARRAY_READERS(int) \ + ARRAY_READERS(unsigned int) \ + ARRAY_READERS(unsigned long int) \ + ARRAY_READERS(double) \ + ARRAY_READERS(float) +#undef ARRAY_READERS + +} // namespace out diff --git a/src/output/utils/readers.h b/src/output/utils/readers.h new file mode 100644 index 000000000..b4fdc8ab0 --- /dev/null +++ b/src/output/utils/readers.h @@ -0,0 +1,48 @@ +/** + * @file output/utils/readers.h + * @brief + * Defines generic reader functions. + * @implements + * - out::ReadVariable<> -> void + * - out::Read1DArray<> -> void + * - out::Read2DArray<> -> void + * @cpp: + * - readers.cpp + * @namespaces: + * - out:: + */ + +#ifndef OUTPUT_UTILS_READERS_H +#define OUTPUT_UTILS_READERS_H + +#include "arch/kokkos_aliases.h" + +#include + +#include + +namespace out { + + template + void ReadVariable(adios2::IO&, adios2::Engine&, const std::string&, T&, std::size_t); + + template + void Read1DArray(adios2::IO&, + adios2::Engine&, + const std::string&, + array_t&, + std::size_t, + std::size_t); + + template + void Read2DArray(adios2::IO&, + adios2::Engine&, + const std::string&, + array_t&, + unsigned short, + std::size_t, + std::size_t); + +} // namespace out + +#endif // OUTPUT_UTILS_READERS_H diff --git a/src/output/utils/writers.cpp b/src/output/utils/writers.cpp new file mode 100644 index 000000000..c7d1bbd74 --- /dev/null +++ b/src/output/utils/writers.cpp @@ -0,0 +1,94 @@ +#include "output/utils/writers.h" + +#include "arch/kokkos_aliases.h" + +#include + +#include + +namespace out { + + template + void WriteVariable(adios2::IO& io, + adios2::Engine& writer, + const std::string& name, + const T& data, + std::size_t global_size, + std::size_t local_offset) { + auto var = io.InquireVariable(name); + var.SetShape({ global_size }); + var.SetSelection(adios2::Box({ local_offset }, { 1 })); + writer.Put(var, &data); + } + + template + void Write1DArray(adios2::IO& io, + adios2::Engine& writer, + const std::string& name, + const array_t& data, + std::size_t local_size, + std::size_t local_offset, + std::size_t global_size) { + const auto slice = range_tuple_t(0, local_size); + auto var = io.InquireVariable(name); + var.SetShape({ global_size }); + var.SetSelection(adios2::Box({ local_offset }, { local_size })); + + auto data_h = Kokkos::create_mirror_view(data); + Kokkos::deep_copy(data_h, data); + auto data_sub = Kokkos::subview(data_h, slice); + writer.Put(var, data_sub.data(), adios2::Mode::Sync); + } + + template + void Write2DArray(adios2::IO& io, + adios2::Engine& writer, + const std::string& name, + const array_t& data, + unsigned short dim2_size, + std::size_t local_size, + std::size_t local_offset, + std::size_t global_size) { + const auto slice = range_tuple_t(0, local_size); + auto var = io.InquireVariable(name); + + var.SetShape({ global_size, dim2_size }); + var.SetSelection( + adios2::Box({ local_offset, 0 }, { local_size, dim2_size })); + + auto data_h = Kokkos::create_mirror_view(data); + Kokkos::deep_copy(data_h, data); + auto data_sub = Kokkos::subview(data_h, slice, range_tuple_t(0, dim2_size)); + writer.Put(var, data_sub.data(), adios2::Mode::Sync); + } + +#define ARRAY_WRITERS(T) \ + template void WriteVariable(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + const T&, \ + std::size_t, \ + std::size_t); \ + template void Write1DArray(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + const array_t&, \ + std::size_t, \ + std::size_t, \ + std::size_t); \ + template void Write2DArray(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + const array_t&, \ + unsigned short, \ + std::size_t, \ + std::size_t, \ + std::size_t); \ + ARRAY_WRITERS(int) \ + ARRAY_WRITERS(unsigned int) \ + ARRAY_WRITERS(unsigned long int) \ + ARRAY_WRITERS(double) \ + ARRAY_WRITERS(float) +#undef ARRAY_WRITERS + +} // namespace out diff --git a/src/output/utils/writers.h b/src/output/utils/writers.h new file mode 100644 index 000000000..3e4bd4cb8 --- /dev/null +++ b/src/output/utils/writers.h @@ -0,0 +1,55 @@ +/** + * @file output/utils/writers.h + * @brief + * Defines generic writer functions. + * @implements + * - out::WriteVariable<> -> void + * - out::Write1DArray<> -> void + * - out::Write2DArray<> -> void + * @cpp: + * - writers.cpp + * @namespaces: + * - out:: + */ + +#ifndef OUTPUT_UTILS_WRITERS_H +#define OUTPUT_UTILS_WRITERS_H + +#include "arch/kokkos_aliases.h" + +#include + +#include + +namespace out { + + template + void WriteVariable(adios2::IO&, + adios2::Engine&, + const std::string&, + const T&, + std::size_t, + std::size_t); + + template + void Write1DArray(adios2::IO&, + adios2::Engine&, + const std::string&, + const array_t&, + std::size_t, + std::size_t, + std::size_t); + + template + void Write2DArray(adios2::IO&, + adios2::Engine&, + const std::string&, + const array_t&, + unsigned short, + std::size_t, + std::size_t, + std::size_t); + +} // namespace out + +#endif // OUTPUT_UTILS_WRITERS_H From 479c0cf49bde984f73bcdd6924b3e945e44cb421 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 21 Oct 2025 16:04:53 -0400 Subject: [PATCH 02/15] hdf5 off by default --- cmake/adios2Config.cmake | 2 +- dev/nix/adios2.nix | 11 ++++++++--- src/global/defaults.h | 2 +- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/cmake/adios2Config.cmake b/cmake/adios2Config.cmake index a4ce46179..0ff3a4e89 100644 --- a/cmake/adios2Config.cmake +++ b/cmake/adios2Config.cmake @@ -13,7 +13,7 @@ set(ADIOS2_USE_Fortran # Format/compression support set(ADIOS2_USE_HDF5 - ON + OFF CACHE BOOL "Use HDF5 for ADIOS2") set(ADIOS2_USE_MPI diff --git a/dev/nix/adios2.nix b/dev/nix/adios2.nix index 0418b71cd..fb574c302 100644 --- a/dev/nix/adios2.nix +++ b/dev/nix/adios2.nix @@ -19,7 +19,7 @@ let BUILD_TESTING = "OFF"; ADIOS2_BUILD_EXAMPLES = "OFF"; ADIOS2_USE_MPI = if mpi then "ON" else "OFF"; - ADIOS2_HAVE_HDF5_VOL = if mpi then "ON" else "OFF"; + ADIOS2_HAVE_HDF5_VOL = if (mpi && hdf5) then "ON" else "OFF"; CMAKE_BUILD_TYPE = "Release"; }; stdenv = pkgs.gcc13Stdenv; @@ -40,8 +40,13 @@ stdenv.mkDerivation { propagatedBuildInputs = [ pkgs.gcc13 - ] ++ (if hdf5 then (if mpi then [ pkgs.hdf5-mpi ] else [ pkgs.hdf5-cpp ]) else [ ]); - # ++ (if mpi then [ pkgs.openmpi ] else [ ]); + ] + ++ ( + if hdf5 then + (if mpi then [ pkgs.hdf5-mpi ] else [ pkgs.hdf5-cpp ]) + else + (if mpi then [ pkgs.mpi ] else [ ]) + ); configurePhase = '' cmake -B build $src ${ diff --git a/src/global/defaults.h b/src/global/defaults.h index 9513493b1..b7ae6547d 100644 --- a/src/global/defaults.h +++ b/src/global/defaults.h @@ -51,7 +51,7 @@ namespace ntt::defaults { } // namespace bc namespace output { - const std::string format = "hdf5"; + const std::string format = "BPFile"; const timestep_t interval = 100; const unsigned short mom_smooth = 0; const npart_t prtl_stride = 100; From c70c3531e213ae33c37d575073f103a4a6ba0458 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 21 Oct 2025 16:07:58 -0400 Subject: [PATCH 03/15] checkpoint functionality for fields placed in fields struct --- src/CMakeLists.txt | 4 - src/checkpoint/CMakeLists.txt | 37 --- src/checkpoint/reader.cpp | 183 ----------- src/checkpoint/reader.h | 57 ---- src/checkpoint/tests/CMakeLists.txt | 30 -- src/checkpoint/tests/checkpoint-mpi.cpp | 273 ----------------- src/checkpoint/tests/checkpoint-nompi.cpp | 248 --------------- src/checkpoint/writer.cpp | 285 ------------------ src/framework/CMakeLists.txt | 7 +- src/framework/containers/fields.h | 16 + src/framework/containers/fields_io.cpp | 96 ++++++ src/framework/containers/particles_io.cpp | 64 ++-- src/framework/domain/checkpoint.cpp | 72 ++--- src/framework/domain/metadomain.h | 2 +- src/framework/parameters.cpp | 22 +- src/framework/parameters.h | 4 + src/output/CMakeLists.txt | 1 + src/output/checkpoint.cpp | 104 +++++++ .../writer.h => output/checkpoint.h} | 30 +- 19 files changed, 323 insertions(+), 1212 deletions(-) delete mode 100644 src/checkpoint/CMakeLists.txt delete mode 100644 src/checkpoint/reader.cpp delete mode 100644 src/checkpoint/reader.h delete mode 100644 src/checkpoint/tests/CMakeLists.txt delete mode 100644 src/checkpoint/tests/checkpoint-mpi.cpp delete mode 100644 src/checkpoint/tests/checkpoint-nompi.cpp delete mode 100644 src/checkpoint/writer.cpp create mode 100644 src/framework/containers/fields_io.cpp create mode 100644 src/output/checkpoint.cpp rename src/{checkpoint/writer.h => output/checkpoint.h} (65%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5f2afc824..31114c330 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,10 +33,6 @@ add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) add_subdirectory(${SRC_DIR}/engines ${CMAKE_CURRENT_BINARY_DIR}/engines) add_subdirectory(${SRC_DIR}/output ${CMAKE_CURRENT_BINARY_DIR}/output) -if(${output}) - add_subdirectory(${SRC_DIR}/checkpoint ${CMAKE_CURRENT_BINARY_DIR}/checkpoint) -endif() - set(ENTITY ${PROJECT_NAME}.xc) set(SOURCES ${SRC_DIR}/entity.cpp) diff --git a/src/checkpoint/CMakeLists.txt b/src/checkpoint/CMakeLists.txt deleted file mode 100644 index 096aad690..000000000 --- a/src/checkpoint/CMakeLists.txt +++ /dev/null @@ -1,37 +0,0 @@ -# cmake-lint: disable=C0103 -# ------------------------------ -# @defines: ntt_checkpoint [STATIC/SHARED] -# -# @sources: -# -# * writer.cpp -# * reader.cpp -# -# @includes: -# -# * ../ -# -# @depends: -# -# * ntt_global [required] -# -# @uses: -# -# * kokkos [required] -# * ADIOS2 [required] -# * mpi [optional] -# ------------------------------ - -set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}) -set(SOURCES ${SRC_DIR}/writer.cpp ${SRC_DIR}/reader.cpp) -add_library(ntt_checkpoint ${SOURCES}) - -set(libs ntt_global) -add_dependencies(ntt_checkpoint ${libs}) -target_link_libraries(ntt_checkpoint PUBLIC ${libs}) -target_link_libraries(ntt_checkpoint PRIVATE stdc++fs) - -target_include_directories( - ntt_checkpoint - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../ - INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/../) diff --git a/src/checkpoint/reader.cpp b/src/checkpoint/reader.cpp deleted file mode 100644 index b5d6f0c44..000000000 --- a/src/checkpoint/reader.cpp +++ /dev/null @@ -1,183 +0,0 @@ -#include "checkpoint/reader.h" - -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "utils/error.h" -#include "utils/formatting.h" -#include "utils/log.h" - -#include -#include - -#if defined(MPI_ENABLED) - #include -#endif - -#include -#include -#include - -namespace checkpoint { - - template - void ReadFields(adios2::IO& io, - adios2::Engine& reader, - const std::string& field, - const adios2::Box& range, - ndfield_t& array) { - logger::Checkpoint(fmt::format("Reading field: %s", field.c_str()), HERE); - auto field_var = io.InquireVariable(field); - if (field_var) { - field_var.SetSelection(range); - - auto array_h = Kokkos::create_mirror_view(array); - reader.Get(field_var, array_h.data(), adios2::Mode::Sync); - Kokkos::deep_copy(array, array_h); - } else { - raise::Error(fmt::format("Field variable: %s not found", field.c_str()), - HERE); - } - } - - auto ReadParticleCount(adios2::IO& io, - adios2::Engine& reader, - spidx_t s, - std::size_t local_dom, - std::size_t ndomains) -> std::pair { - logger::Checkpoint(fmt::format("Reading particle count for: %d", s + 1), HERE); - auto npart_var = io.InquireVariable(fmt::format("s%d_npart", s + 1)); - if (npart_var) { - raise::ErrorIf(npart_var.Shape()[0] != ndomains, - "npart_var.Shape()[0] != ndomains", - HERE); - raise::ErrorIf(npart_var.Shape().size() != 1, - "npart_var.Shape().size() != 1", - HERE); - npart_var.SetSelection(adios2::Box({ local_dom }, { 1 })); - npart_t npart; - reader.Get(npart_var, &npart, adios2::Mode::Sync); - const auto loc_npart = npart; -#if !defined(MPI_ENABLED) - npart_t offset_npart = 0; -#else - std::vector glob_nparts(ndomains); - MPI_Allgather(&loc_npart, - 1, - mpi::get_type(), - glob_nparts.data(), - 1, - mpi::get_type(), - MPI_COMM_WORLD); - npart_t offset_npart = 0; - for (auto d { 0u }; d < local_dom; ++d) { - offset_npart += glob_nparts[d]; - } -#endif - return { loc_npart, offset_npart }; - } else { - raise::Error("npart_var is not found", HERE); - return { 0, 0 }; - } - } - - template - void ReadParticleData(adios2::IO& io, - adios2::Engine& reader, - const std::string& quantity, - spidx_t s, - array_t& array, - npart_t count, - npart_t offset) { - logger::Checkpoint( - fmt::format("Reading quantity: s%d_%s", s + 1, quantity.c_str()), - HERE); - auto var = io.InquireVariable( - fmt::format("s%d_%s", s + 1, quantity.c_str())); - if (var) { - var.SetSelection(adios2::Box({ offset }, { count })); - const auto slice = range_tuple_t(0, count); - auto array_h = Kokkos::create_mirror_view(array); - reader.Get(var, Kokkos::subview(array_h, slice).data(), adios2::Mode::Sync); - Kokkos::deep_copy(Kokkos::subview(array, slice), - Kokkos::subview(array_h, slice)); - } else { - raise::Error( - fmt::format("Variable: s%d_%s not found", s + 1, quantity.c_str()), - HERE); - } - } - - template - void ReadParticlePayloads(adios2::IO& io, - adios2::Engine& reader, - const std::string& suffix, - spidx_t s, - array_t& array, - std::size_t nplds, - npart_t count, - npart_t offset) { - logger::Checkpoint( - fmt::format("Reading quantity: s%d_pld_%s", s + 1, suffix.c_str()), - HERE); - auto var = io.InquireVariable( - fmt::format("s%d_pld_%s", s + 1, suffix.c_str())); - if (var) { - var.SetSelection(adios2::Box({ offset, 0 }, { count, nplds })); - const auto slice = range_tuple_t(0, count); - auto array_h = Kokkos::create_mirror_view(array); - reader.Get(var, - Kokkos::subview(array_h, slice, range_tuple_t(0, nplds)).data(), - adios2::Mode::Sync); - Kokkos::deep_copy(array, array_h); - } else { - raise::Error( - fmt::format("Variable: s%d_pld_%s not found", s + 1, suffix.c_str()), - HERE); - } - } - -#define CHECKPOINT_FIELDS(D, N) \ - template void ReadFields(adios2::IO&, \ - adios2::Engine&, \ - const std::string&, \ - const adios2::Box&, \ - ndfield_t&); - CHECKPOINT_FIELDS(Dim::_1D, 3) - CHECKPOINT_FIELDS(Dim::_2D, 3) - CHECKPOINT_FIELDS(Dim::_3D, 3) - CHECKPOINT_FIELDS(Dim::_1D, 6) - CHECKPOINT_FIELDS(Dim::_2D, 6) - CHECKPOINT_FIELDS(Dim::_3D, 6) -#undef CHECKPOINT_FIELDS - -#define CHECKPOINT_PARTICLE_DATA(T) \ - template void ReadParticleData(adios2::IO&, \ - adios2::Engine&, \ - const std::string&, \ - spidx_t, \ - array_t&, \ - npart_t, \ - npart_t); - CHECKPOINT_PARTICLE_DATA(int) - CHECKPOINT_PARTICLE_DATA(float) - CHECKPOINT_PARTICLE_DATA(double) - CHECKPOINT_PARTICLE_DATA(short) -#undef CHECKPOINT_PARTICLE_DATA - -#define CHECKPOINT_PARTICLE_PAYLOADS(T) \ - template void ReadParticlePayloads(adios2::IO&, \ - adios2::Engine&, \ - const std::string&, \ - spidx_t, \ - array_t&, \ - std::size_t, \ - npart_t, \ - npart_t); - CHECKPOINT_PARTICLE_PAYLOADS(int) - CHECKPOINT_PARTICLE_PAYLOADS(float) - CHECKPOINT_PARTICLE_PAYLOADS(double) - CHECKPOINT_PARTICLE_PAYLOADS(npart_t) -#undef CHECKPOINT_PARTICLE_DATA - -} // namespace checkpoint diff --git a/src/checkpoint/reader.h b/src/checkpoint/reader.h deleted file mode 100644 index 78cc73eb7..000000000 --- a/src/checkpoint/reader.h +++ /dev/null @@ -1,57 +0,0 @@ -/** - * @file checkpoint/reader.h - * @brief Function for reading field & particle data from checkpoint files - * @implements - * - checkpoint::ReadFields -> void - * - checkpoint::ReadParticleData -> void - * - checkpoint::ReadParticleCount -> std::pair - * @cpp: - * - reader.cpp - * @namespaces: - * - checkpoint:: - */ - -#ifndef CHECKPOINT_READER_H -#define CHECKPOINT_READER_H - -#include "arch/kokkos_aliases.h" - -#include - -#include -#include - -namespace checkpoint { - - template - void ReadFields(adios2::IO&, - adios2::Engine&, - const std::string&, - const adios2::Box&, - ndfield_t&); - - auto ReadParticleCount(adios2::IO&, adios2::Engine&, spidx_t, std::size_t, std::size_t) - -> std::pair; - - template - void ReadParticleData(adios2::IO&, - adios2::Engine&, - const std::string&, - spidx_t, - array_t&, - npart_t, - npart_t); - - template - void ReadParticlePayloads(adios2::IO&, - adios2::Engine&, - const std::string&, - spidx_t, - array_t&, - std::size_t, - npart_t, - npart_t); - -} // namespace checkpoint - -#endif // CHECKPOINT_READER_H diff --git a/src/checkpoint/tests/CMakeLists.txt b/src/checkpoint/tests/CMakeLists.txt deleted file mode 100644 index cbfd63aa9..000000000 --- a/src/checkpoint/tests/CMakeLists.txt +++ /dev/null @@ -1,30 +0,0 @@ -# cmake-lint: disable=C0103,C0111 -# ------------------------------ -# @brief: Generates tests for the `ntt_checkpoint` module -# -# @uses: -# -# * kokkos [required] -# * adios2 [required] -# * mpi [optional] -# ------------------------------ - -set(SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../) - -function(gen_test title) - set(exec test-output-${title}.xc) - set(src ${title}.cpp) - add_executable(${exec} ${src}) - - set(libs ntt_checkpoint ntt_global) - add_dependencies(${exec} ${libs}) - target_link_libraries(${exec} PRIVATE ${libs} stdc++fs) - - add_test(NAME "CHECKPOINT::${title}" COMMAND "${exec}") -endfunction() - -if(NOT ${mpi}) - gen_test(checkpoint-nompi) -else() - gen_test(checkpoint-mpi) -endif() diff --git a/src/checkpoint/tests/checkpoint-mpi.cpp b/src/checkpoint/tests/checkpoint-mpi.cpp deleted file mode 100644 index 2372d81bc..000000000 --- a/src/checkpoint/tests/checkpoint-mpi.cpp +++ /dev/null @@ -1,273 +0,0 @@ -#include "enums.h" -#include "global.h" - -#include "utils/comparators.h" - -#include "checkpoint/reader.h" -#include "checkpoint/writer.h" - -#include -#include -#include -#include - -#include -#include -#include - -using namespace ntt; -using namespace checkpoint; - -void cleanup() { - namespace fs = std::filesystem; - fs::path temp_path { "chck" }; - fs::remove_all(temp_path); -} - -auto main(int argc, char* argv[]) -> int { - Kokkos::initialize(argc, argv); - MPI_Init(&argc, &argv); - int rank, size; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &size); - - try { - // assuming 4 ranks - // |------|------| - // | 2 | 3 | - // |------|------| - // | | | - // | 0 | 1 | - // |------|------| - const std::size_t g_nx1 = 20; - const std::size_t g_nx2 = 15; - const std::size_t g_nx1_gh = g_nx1 + 4 * N_GHOSTS; - const std::size_t g_nx2_gh = g_nx2 + 4 * N_GHOSTS; - - const std::size_t l_nx1 = 10; - const std::size_t l_nx2 = (rank < 2) ? 10 : 5; - - const std::size_t l_nx1_gh = l_nx1 + 2 * N_GHOSTS; - const std::size_t l_nx2_gh = l_nx2 + 2 * N_GHOSTS; - - const std::size_t l_corner_x1 = (rank % 2 == 0) ? 0 : l_nx1_gh; - const std::size_t l_corner_x2 = (rank < 2) ? 0 : l_nx2_gh; - - const std::size_t i1min = N_GHOSTS; - const std::size_t i2min = N_GHOSTS; - const std::size_t i1max = l_nx1 + N_GHOSTS; - const std::size_t i2max = l_nx2 + N_GHOSTS; - - const std::size_t npart1 = (rank % 2 + rank) * 23 + 100; - const std::size_t npart2 = (rank % 2 + rank) * 37 + 100; - - std::size_t npart1_offset = 0; - std::size_t npart2_offset = 0; - - std::size_t npart1_globtot = 0; - std::size_t npart2_globtot = 0; - - for (auto r = 0; r < rank - 1; ++r) { - npart1_offset += (r % 2 + r) * 23 + 100; - npart2_offset += (r % 2 + r) * 37 + 100; - } - - for (auto r = 0; r < size; ++r) { - npart1_globtot += (r % 2 + r) * 23 + 100; - npart2_globtot += (r % 2 + r) * 37 + 100; - } - - // init data - ndfield_t field1 { "fld1", l_nx1_gh, l_nx2_gh }; - ndfield_t field2 { "fld2", l_nx1_gh, l_nx2_gh }; - - array_t i1 { "i_1", npart1 }; - array_t u1 { "u_1", npart1 }; - array_t i2 { "i_2", npart2 }; - array_t u2 { "u_2", npart2 }; - array_t plds1 { "plds_1", npart1, 3 }; - - { - // fill data - Kokkos::parallel_for( - "fillFlds", - CreateRangePolicy({ i1min, i2min }, { i1max, i2max }), - Lambda(index_t i1, index_t i2) { - field1(i1, i2, 0) = static_cast(i1 + i2); - field1(i1, i2, 1) = static_cast(i1 * i2); - field1(i1, i2, 2) = static_cast(i1 / i2); - field1(i1, i2, 3) = static_cast(i1 - i2); - field1(i1, i2, 4) = static_cast(i2 / i1); - field1(i1, i2, 5) = static_cast(i1); - field2(i1, i2, 0) = static_cast(-(i1 + i2)); - field2(i1, i2, 1) = static_cast(-(i1 * i2)); - field2(i1, i2, 2) = static_cast(-(i1 / i2)); - field2(i1, i2, 3) = static_cast(-(i1 - i2)); - field2(i1, i2, 4) = static_cast(-(i2 / i1)); - field2(i1, i2, 5) = static_cast(-i1); - }); - Kokkos::parallel_for( - "fillPrtl1", - npart1, - Lambda(index_t p) { - u1(p) = static_cast(p); - i1(p) = static_cast(p); - plds1(p, 0) = static_cast(p); - plds1(p, 1) = static_cast(p * p); - plds1(p, 2) = static_cast(p * p * p); - }); - Kokkos::parallel_for( - "fillPrtl2", - npart2, - Lambda(index_t p) { - u2(p) = -static_cast(p); - i2(p) = -static_cast(p); - }); - } - - adios2::ADIOS adios; - const path_t checkpoint_path { "chck" }; - - { - // write checkpoint - Writer writer; - writer.init(&adios, checkpoint_path, 0, 0.0, 1); - - writer.defineFieldVariables(SimEngine::GRPIC, - { g_nx1_gh, g_nx2_gh }, - { l_corner_x1, l_corner_x2 }, - { l_nx1_gh, l_nx2_gh }); - - writer.defineParticleVariables(Coord::Sph, Dim::_2D, 2, { 3, 0 }); - - writer.beginSaving(0, 0.0); - - writer.saveField("em", field1); - writer.saveField("em0", field2); - - writer.savePerDomainVariable("s1_npart", 1, 0, npart1); - writer.savePerDomainVariable("s2_npart", 1, 0, npart2); - - writer.saveParticleQuantity("s1_i1", - npart1_globtot, - npart1_offset, - npart1, - i1); - writer.saveParticleQuantity("s1_ux1", - npart1_globtot, - npart1_offset, - npart1, - u1); - writer.saveParticleQuantity("s2_i1", - npart2_globtot, - npart2_offset, - npart2, - i2); - writer.saveParticleQuantity("s2_ux1", - npart2_globtot, - npart2_offset, - npart2, - u2); - - writer.saveParticlePayloads("s1_plds", - 3, - npart1_globtot, - npart1_offset, - npart1, - plds1); - - writer.endSaving(); - } - - { - // read checkpoint - ndfield_t field1_read { "fld1_read", l_nx1_gh, l_nx2_gh }; - ndfield_t field2_read { "fld2_read", l_nx1_gh, l_nx2_gh }; - - array_t i1_read { "i_1", npart1 }; - array_t u1_read { "u_1", npart1 }; - array_t i2_read { "i_2", npart2 }; - array_t u2_read { "u_2", npart2 }; - array_t plds1_read { "plds_1", npart1, 3 }; - - adios2::IO io = adios.DeclareIO("checkpointRead"); - adios2::Engine reader = io.Open(checkpoint_path / "step-00000000.bp", - adios2::Mode::Read); - reader.BeginStep(); - - auto fieldRange = adios2::Box({ l_corner_x1, l_corner_x2, 0 }, - { l_nx1_gh, l_nx2_gh, 6 }); - ReadFields(io, reader, "em", fieldRange, field1_read); - ReadFields(io, reader, "em0", fieldRange, field2_read); - - auto [nprtl1, noff1] = ReadParticleCount(io, reader, 0, rank, size); - auto [nprtl2, noff2] = ReadParticleCount(io, reader, 1, rank, size); - - ReadParticleData(io, reader, "ux1", 0, u1_read, nprtl1, noff1); - ReadParticleData(io, reader, "ux1", 1, u2_read, nprtl2, noff2); - ReadParticleData(io, reader, "i1", 0, i1_read, nprtl1, noff1); - ReadParticleData(io, reader, "i1", 1, i2_read, nprtl2, noff2); - ReadParticlePayloads(io, reader, 0, plds1_read, 3, nprtl1, noff1); - - reader.EndStep(); - reader.Close(); - - // check the validity - Kokkos::parallel_for( - "checkFields", - CreateRangePolicy({ 0, 0 }, { l_nx1_gh, l_nx2_gh }), - Lambda(index_t i1, index_t i2) { - for (int i = 0; i < 6; ++i) { - if (not cmp::AlmostEqual(field1(i1, i2, i), field1_read(i1, i2, i))) { - raise::KernelError(HERE, "Field1 read failed"); - } - if (not cmp::AlmostEqual(field2(i1, i2, i), field2_read(i1, i2, i))) { - raise::KernelError(HERE, "Field2 read failed"); - } - } - }); - - raise::ErrorIf(npart1 != nprtl1, "Particle count 1 mismatch", HERE); - raise::ErrorIf(npart2 != nprtl2, "Particle count 2 mismatch", HERE); - raise::ErrorIf(noff1 != npart1_offset, "Particle offset 1 mismatch", HERE); - raise::ErrorIf(noff2 != npart2_offset, "Particle offset 2 mismatch", HERE); - - Kokkos::parallel_for( - "checkPrtl1", - nprtl1, - Lambda(index_t p) { - if (not cmp::AlmostEqual(u1(p), u1_read(p))) { - raise::KernelError(HERE, "u1 read failed"); - } - if (i1(p) != i1_read(p)) { - raise::KernelError(HERE, "i1 read failed"); - } - for (auto l = 0; l < 3; ++l) { - if (not cmp::AlmostEqual(plds1(p, l), plds1_read(p, l))) { - raise::KernelError(HERE, "plds1 read failed"); - } - } - }); - Kokkos::parallel_for( - "checkPrtl2", - nprtl2, - Lambda(index_t p) { - if (not cmp::AlmostEqual(u2(p), u2_read(p))) { - raise::KernelError(HERE, "u2 read failed"); - } - if (i2(p) != i2_read(p)) { - raise::KernelError(HERE, "i2 read failed"); - } - }); - } - - } catch (std::exception& e) { - std::cerr << e.what() << std::endl; - cleanup(); - Kokkos::finalize(); - return 1; - } - cleanup(); - Kokkos::finalize(); - return 0; -} diff --git a/src/checkpoint/tests/checkpoint-nompi.cpp b/src/checkpoint/tests/checkpoint-nompi.cpp deleted file mode 100644 index 7be41ca3f..000000000 --- a/src/checkpoint/tests/checkpoint-nompi.cpp +++ /dev/null @@ -1,248 +0,0 @@ -#include "enums.h" -#include "global.h" - -#include "utils/comparators.h" - -#include "checkpoint/reader.h" -#include "checkpoint/writer.h" - -#include -#include -#include - -#include -#include - -using namespace ntt; -using namespace checkpoint; - -void cleanup() { - namespace fs = std::filesystem; - fs::path temp_path { "chck" }; - fs::remove_all(temp_path); -} - -auto main(int argc, char* argv[]) -> int { - Kokkos::initialize(argc, argv); - - try { - constexpr auto nx1 = 10; - constexpr auto nx1_gh = nx1 + 2 * N_GHOSTS; - constexpr auto nx2 = 13; - constexpr auto nx2_gh = nx2 + 2 * N_GHOSTS; - constexpr auto nx3 = 9; - constexpr auto nx3_gh = nx3 + 2 * N_GHOSTS; - constexpr auto i1min = N_GHOSTS; - constexpr auto i2min = N_GHOSTS; - constexpr auto i3min = N_GHOSTS; - constexpr auto i1max = nx1 + N_GHOSTS; - constexpr auto i2max = nx2 + N_GHOSTS; - constexpr auto i3max = nx3 + N_GHOSTS; - constexpr auto npart1 = 100; - constexpr auto npart2 = 100; - - // init data - ndfield_t field1 { "fld1", nx1_gh, nx2_gh, nx3_gh }; - ndfield_t field2 { "fld2", nx1_gh, nx2_gh, nx3_gh }; - - array_t i1 { "i_1", npart1 }; - array_t u1 { "u_1", npart1 }; - array_t i2 { "i_2", npart2 }; - array_t u2 { "u_2", npart2 }; - - array_t pldr_2 { "pldr_2", npart2, 2 }; - - array_t pldi_1 { "pldi_1", npart1, 1 }; - array_t pldi_2 { "pldi_2", npart2, 2 }; - - { - // fill data - Kokkos::parallel_for( - "fillFlds", - CreateRangePolicy({ i1min, i2min, i3min }, - { i1max, i2max, i3max }), - Lambda(index_t i1, index_t i2, index_t i3) { - const auto i1_ = static_cast(i1); - const auto i2_ = static_cast(i2); - const auto i3_ = static_cast(i3); - field1(i1, i2, i3, 0) = i1_ + i2_ + i3_; - field1(i1, i2, i3, 1) = i1_ * i2_ / i3_; - field1(i1, i2, i3, 2) = i1_ / i2_ * i3_; - field1(i1, i2, i3, 3) = i1_ + i2_ - i3_; - field1(i1, i2, i3, 4) = i1_ * i2_ + i3_; - field1(i1, i2, i3, 5) = i1_ / i2_ - i3_; - field2(i1, i2, i3, 0) = -(i1_ + i2_ + i3_); - field2(i1, i2, i3, 1) = -(i1_ * i2_ / i3_); - field2(i1, i2, i3, 2) = -(i1_ / i2_ * i3_); - field2(i1, i2, i3, 3) = -(i1_ + i2_ - i3_); - field2(i1, i2, i3, 4) = -(i1_ * i2_ + i3_); - field2(i1, i2, i3, 5) = -(i1_ / i2_ - i3_); - }); - Kokkos::parallel_for( - "fillPrtl1", - npart1, - Lambda(index_t p) { - u1(p) = static_cast(p); - i1(p) = static_cast(p); - pldi_1(p, 0) = static_cast(p * 10); - }); - Kokkos::parallel_for( - "fillPrtl2", - npart2, - Lambda(index_t p) { - u2(p) = -static_cast(p); - i2(p) = -static_cast(p); - pldr_2(p, 0) = static_cast(p); - pldr_2(p, 1) = static_cast(p * 2); - pldi_2(p, 0) = static_cast(p * 3); - pldi_2(p, 1) = static_cast(p * 4); - }); - } - - adios2::ADIOS adios; - const path_t checkpoint_path { "chck" }; - - { - // write checkpoint - Writer writer {}; - writer.init(&adios, checkpoint_path, 0, 0.0, 1); - - writer.defineFieldVariables(SimEngine::GRPIC, - { nx1_gh, nx2_gh, nx3_gh }, - { 0, 0, 0 }, - { nx1_gh, nx2_gh, nx3_gh }); - writer.defineParticleVariables(Coord::Sph, Dim::_3D, 2, { 0, 2 }, { 1, 2 }); - - writer.beginSaving(0, 0.0); - - writer.saveField("em", field1); - writer.saveField("em0", field2); - - writer.savePerDomainVariable("s1_npart", 1, 0, npart1); - writer.savePerDomainVariable("s2_npart", 1, 0, npart2); - - writer.saveParticleQuantity("s1_i1", npart1, 0, npart1, i1); - writer.saveParticleQuantity("s1_ux1", npart1, 0, npart1, u1); - writer.saveParticleQuantity("s2_i1", npart2, 0, npart2, i2); - writer.saveParticleQuantity("s2_ux1", npart2, 0, npart2, u2); - - writer.saveParticlePayloads("s2_pld_r", 2, npart2, 0, npart2, pldr_2); - - writer.saveParticlePayloads("s1_pld_i", 1, npart1, 0, npart1, pldi_1); - writer.saveParticlePayloads("s2_pld_i", 2, npart2, 0, npart2, pldi_2); - - writer.endSaving(); - } - - { - // read checkpoint - ndfield_t field1_read { "fld1_read", nx1_gh, nx2_gh, nx3_gh }; - ndfield_t field2_read { "fld2_read", nx1_gh, nx2_gh, nx3_gh }; - - array_t i1_read { "i_1", npart1 }; - array_t u1_read { "u_1", npart1 }; - array_t i2_read { "i_2", npart2 }; - array_t u2_read { "u_2", npart2 }; - - array_t pldr_2_read { "pldr_2", npart2, 2 }; - - array_t pldi_1_read { "pldi_1", npart1, 1 }; - array_t pldi_2_read { "pldi_2", npart2, 2 }; - - adios2::IO io = adios.DeclareIO("checkpointRead"); - adios2::Engine reader = io.Open(checkpoint_path / "step-00000000.bp", - adios2::Mode::Read); - reader.BeginStep(); - - auto fieldRange = adios2::Box({ 0, 0, 0, 0 }, - { nx1_gh, nx2_gh, nx3_gh, 6 }); - ReadFields(io, reader, "em", fieldRange, field1_read); - ReadFields(io, reader, "em0", fieldRange, field2_read); - - auto [nprtl1, noff1] = ReadParticleCount(io, reader, 0, 0, 1); - auto [nprtl2, noff2] = ReadParticleCount(io, reader, 1, 0, 1); - - ReadParticleData(io, reader, "ux1", 0, u1_read, nprtl1, noff1); - ReadParticleData(io, reader, "ux1", 1, u2_read, nprtl2, noff2); - ReadParticleData(io, reader, "i1", 0, i1_read, nprtl1, noff1); - ReadParticleData(io, reader, "i1", 1, i2_read, nprtl2, noff2); - - ReadParticlePayloads(io, reader, "r", 1, pldr_2_read, 2, nprtl2, noff2); - - ReadParticlePayloads(io, reader, "i", 0, pldi_1_read, 1, nprtl1, noff1); - ReadParticlePayloads(io, reader, "i", 1, pldi_2_read, 2, nprtl2, noff2); - - reader.EndStep(); - reader.Close(); - - // check the validity - Kokkos::parallel_for( - "checkFields", - CreateRangePolicy({ 0, 0, 0 }, { nx1_gh, nx2_gh, nx3_gh }), - Lambda(index_t i1, index_t i2, index_t i3) { - for (int i = 0; i < 6; ++i) { - if (not cmp::AlmostEqual(field1(i1, i2, i3, i), - field1_read(i1, i2, i3, i))) { - raise::KernelError(HERE, "Field1 read failed"); - } - if (not cmp::AlmostEqual(field2(i1, i2, i3, i), - field2_read(i1, i2, i3, i))) { - raise::KernelError(HERE, "Field2 read failed"); - } - } - }); - - raise::ErrorIf(npart1 != nprtl1, "Particle count 1 mismatch", HERE); - raise::ErrorIf(npart2 != nprtl2, "Particle count 2 mismatch", HERE); - raise::ErrorIf(noff1 != 0, "Particle offset 1 mismatch", HERE); - raise::ErrorIf(noff2 != 0, "Particle offset 2 mismatch", HERE); - - Kokkos::parallel_for( - "checkPrtl1", - npart1, - Lambda(index_t p) { - if (not cmp::AlmostEqual(u1(p), u1_read(p))) { - raise::KernelError(HERE, "u1 read failed"); - } - if (i1(p) != i1_read(p)) { - raise::KernelError(HERE, "i1 read failed"); - } - if (pldi_1(p, 0) != pldi_1_read(p, 0)) { - raise::KernelError(HERE, "pldi_1 read failed"); - } - }); - Kokkos::parallel_for( - "checkPrtl2", - npart2, - Lambda(index_t p) { - if (not cmp::AlmostEqual(u2(p), u2_read(p))) { - raise::KernelError(HERE, "u2 read failed"); - } - if (i2(p) != i2_read(p)) { - raise::KernelError(HERE, "i2 read failed"); - } - if (not cmp::AlmostEqual(pldr_2(p, 0), pldr_2_read(p, 0))) { - raise::KernelError(HERE, "pldr_2(0) read failed"); - } - if (not cmp::AlmostEqual(pldr_2(p, 1), pldr_2_read(p, 1))) { - raise::KernelError(HERE, "pldr_2(1) read failed"); - } - if (pldi_2(p, 0) != pldi_2_read(p, 0)) { - raise::KernelError(HERE, "pldi_2(0) read failed"); - } - if (pldi_2(p, 1) != pldi_2_read(p, 1)) { - raise::KernelError(HERE, "pldi_2(1) read failed"); - } - }); - } - - } catch (std::exception& e) { - std::cerr << e.what() << std::endl; - cleanup(); - Kokkos::finalize(); - return 1; - } - cleanup(); - Kokkos::finalize(); - return 0; -} diff --git a/src/checkpoint/writer.cpp b/src/checkpoint/writer.cpp deleted file mode 100644 index ebb8663ae..000000000 --- a/src/checkpoint/writer.cpp +++ /dev/null @@ -1,285 +0,0 @@ -#include "checkpoint/writer.h" - -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "utils/error.h" -#include "utils/formatting.h" -#include "utils/log.h" - -#include "framework/parameters.h" - -#include -#include - -#include -#include -#include -#include - -namespace checkpoint { - - void Writer::init(adios2::ADIOS* ptr_adios, - const path_t& checkpoint_root, - timestep_t interval, - simtime_t interval_time, - int keep, - const std::string& walltime) { - m_keep = keep; - m_checkpoint_root = checkpoint_root; - m_enabled = keep != 0; - if (not m_enabled) { - return; - } - m_tracker.init("checkpoint", interval, interval_time, walltime); - p_adios = ptr_adios; - raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE); - - m_io = p_adios->DeclareIO("Entity::Checkpoint"); - m_io.SetEngine("BPFile"); - - m_io.DefineVariable("Step"); - m_io.DefineVariable("Time"); - m_io.DefineAttribute("NGhosts", ntt::N_GHOSTS); - - CallOnce( - [](auto&& checkpoint_root) { - if (!std::filesystem::exists(checkpoint_root)) { - std::filesystem::create_directory(checkpoint_root); - } - }, - m_checkpoint_root); - } - - void Writer::defineFieldVariables(const ntt::SimEngine& S, - const std::vector& glob_shape, - const std::vector& loc_corner, - const std::vector& loc_shape) { - auto gs6 = std::vector(glob_shape.begin(), glob_shape.end()); - auto lc6 = std::vector(loc_corner.begin(), loc_corner.end()); - auto ls6 = std::vector(loc_shape.begin(), loc_shape.end()); - gs6.push_back(6); - lc6.push_back(0); - ls6.push_back(6); - - m_io.DefineVariable("em", gs6, lc6, ls6); - if (S == ntt::SimEngine::GRPIC) { - m_io.DefineVariable("em0", gs6, lc6, ls6); - auto gs3 = std::vector(glob_shape.begin(), glob_shape.end()); - auto lc3 = std::vector(loc_corner.begin(), loc_corner.end()); - auto ls3 = std::vector(loc_shape.begin(), loc_shape.end()); - gs3.push_back(3); - lc3.push_back(0); - ls3.push_back(3); - m_io.DefineVariable("cur0", gs3, lc3, ls3); - } - } - - // void Writer::defineParticleVariables(const ntt::Coord& C, - // Dimension dim, - // std::size_t nspec, - // const std::vector& npld_r, - // const std::vector& npld_i) { - // raise::ErrorIf( - // npld_r.size() != nspec, - // "Number of real payloads does not match the number of species", - // HERE); - // raise::ErrorIf( - // npld_i.size() != nspec, - // "Number of int payloads does not match the number of species", - // HERE); - // for (auto s { 0u }; s < nspec; ++s) { - // m_io.DefineVariable(fmt::format("s%d_npart", s + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // m_io.DefineVariable(fmt::format("s%d_counter", s + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // - // for (auto d { 0u }; d < dim; ++d) { - // m_io.DefineVariable(fmt::format("s%d_i%d", s + 1, d + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // m_io.DefineVariable(fmt::format("s%d_dx%d", s + 1, d + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // m_io.DefineVariable(fmt::format("s%d_i%d_prev", s + 1, d + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // m_io.DefineVariable(fmt::format("s%d_dx%d_prev", s + 1, d + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // } - // - // if (dim == Dim::_2D and C != ntt::Coord::Cart) { - // m_io.DefineVariable(fmt::format("s%d_phi", s + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // } - // - // for (auto d { 0u }; d < 3; ++d) { - // m_io.DefineVariable(fmt::format("s%d_ux%d", s + 1, d + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // } - // - // m_io.DefineVariable(fmt::format("s%d_tag", s + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // m_io.DefineVariable(fmt::format("s%d_weight", s + 1), - // { adios2::UnknownDim }, - // { adios2::UnknownDim }, - // { adios2::UnknownDim }); - // if (npld_r[s] > 0) { - // m_io.DefineVariable(fmt::format("s%d_pld_r", s + 1), - // { adios2::UnknownDim, npld_r[s] }, - // { adios2::UnknownDim, 0 }, - // { adios2::UnknownDim, npld_r[s] }); - // } - // if (npld_i[s] > 0) { - // m_io.DefineVariable(fmt::format("s%d_pld_i", s + 1), - // { adios2::UnknownDim, npld_i[s] }, - // { adios2::UnknownDim, 0 }, - // { adios2::UnknownDim, npld_i[s] }); - // } - // } - // } - - auto Writer::shouldSave(timestep_t step, simtime_t time) -> bool { - return m_enabled and m_tracker.shouldWrite(step, time); - } - - void Writer::beginSaving(timestep_t step, simtime_t time) { - raise::ErrorIf(!m_enabled, "Checkpoint is not enabled", HERE); - raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE); - if (m_writing_mode) { - raise::Fatal("Already writing", HERE); - } - m_writing_mode = true; - try { - const auto filename = m_checkpoint_root / fmt::format("step-%08lu.bp", step); - const auto metafilename = m_checkpoint_root / - fmt::format("meta-%08lu.toml", step); - m_writer = m_io.Open(filename, adios2::Mode::Write); - m_written.push_back({ filename, metafilename }); - logger::Checkpoint(fmt::format("Writing checkpoint to %s and %s", - filename.c_str(), - metafilename.c_str()), - HERE); - } catch (std::exception& e) { - raise::Fatal(e.what(), HERE); - } - - m_writer.BeginStep(); - m_writer.Put(m_io.InquireVariable("Step"), &step); - m_writer.Put(m_io.InquireVariable("Time"), &time); - } - - void Writer::endSaving() { - raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE); - if (!m_writing_mode) { - raise::Fatal("Not writing", HERE); - } - m_writing_mode = false; - m_writer.EndStep(); - m_writer.Close(); - - // optionally remove the oldest checkpoint - CallOnce([&]() { - if (m_keep > 0 and m_written.size() > (std::size_t)m_keep) { - const auto oldest = m_written.front(); - if (std::filesystem::exists(oldest.first) and - std::filesystem::exists(oldest.second)) { - std::filesystem::remove_all(oldest.first); - std::filesystem::remove(oldest.second); - m_written.erase(m_written.begin()); - } else { - raise::Warning("Checkpoint file does not exist for some reason", HERE); - } - } - }); - } - - void Writer::saveAttrs(const ntt::SimulationParams& params, simtime_t time) { - CallOnce([&]() { - std::ofstream metadata; - if (m_written.empty()) { - raise::Fatal("No checkpoint file to save metadata", HERE); - } - metadata.open(m_written.back().second.c_str()); - metadata << "[metadata]\n" - << " time = " << time << "\n\n" - << params.data() << std::endl; - metadata.close(); - }); - } - - template - void Writer::saveField(const std::string& fieldname, - const ndfield_t& field) { - auto field_h = Kokkos::create_mirror_view(field); - Kokkos::deep_copy(field_h, field); - m_writer.Put(m_io.InquireVariable(fieldname), - field_h.data(), - adios2::Mode::Sync); - } - - // template - // void Writer::saveParticleQuantity(const std::string& quantity, - // npart_t glob_total, - // npart_t loc_offset, - // npart_t loc_size, - // const array_t& data) { - // const auto slice = range_tuple_t(0, loc_size); - // auto var = m_io.InquireVariable(quantity); - // - // var.SetShape({ glob_total }); - // var.SetSelection(adios2::Box({ loc_offset }, { loc_size })); - // - // auto data_h = Kokkos::create_mirror_view(data); - // Kokkos::deep_copy(data_h, data); - // auto data_sub = Kokkos::subview(data_h, slice); - // m_writer.Put(var, data_sub.data(), adios2::Mode::Sync); - // } - // - // template - // void Writer::saveParticlePayloads(const std::string& quantity, - // std::size_t nplds, - // npart_t glob_total, - // npart_t loc_offset, - // npart_t loc_size, - // const array_t& data) { - // const auto slice = range_tuple_t(0, loc_size); - // auto var = m_io.InquireVariable(quantity); - // - // var.SetShape({ glob_total, nplds }); - // var.SetSelection( - // adios2::Box({ loc_offset, 0 }, { loc_size, nplds })); - // - // auto data_h = Kokkos::create_mirror_view(data); - // Kokkos::deep_copy(data_h, data); - // auto data_sub = Kokkos::subview(data_h, slice, range_tuple_t(0, nplds)); - // m_writer.Put(var, data_sub.data(), adios2::Mode::Sync); - // } - -#define CHECKPOINT_FIELD(D, N) \ - template void Writer::saveField(const std::string&, \ - const ndfield_t&); - CHECKPOINT_FIELD(Dim::_1D, 3) - CHECKPOINT_FIELD(Dim::_1D, 6) - CHECKPOINT_FIELD(Dim::_2D, 3) - CHECKPOINT_FIELD(Dim::_2D, 6) - CHECKPOINT_FIELD(Dim::_3D, 3) - CHECKPOINT_FIELD(Dim::_3D, 6) -#undef CHECKPOINT_FIELD - -} // namespace checkpoint diff --git a/src/framework/CMakeLists.txt b/src/framework/CMakeLists.txt index 07328fff6..014967870 100644 --- a/src/framework/CMakeLists.txt +++ b/src/framework/CMakeLists.txt @@ -47,14 +47,15 @@ set(SOURCES if(${output}) list(APPEND SOURCES ${SRC_DIR}/domain/output.cpp) list(APPEND SOURCES ${SRC_DIR}/domain/checkpoint.cpp) + list(APPEND SOURCES ${SRC_DIR}/containers/fields_io.cpp) list(APPEND SOURCES ${SRC_DIR}/containers/particles_io.cpp) endif() +if(${mpi}) + list(APPEND SOURCES ${SRC_DIR}/containers/particles_comm.cpp) +endif() add_library(ntt_framework ${SOURCES}) set(libs ntt_global ntt_metrics ntt_kernels ntt_output) -if(${output}) - list(APPEND libs ntt_checkpoint) -endif() add_dependencies(ntt_framework ${libs}) target_link_libraries(ntt_framework PUBLIC ${libs}) target_link_libraries(ntt_framework PRIVATE stdc++fs) diff --git a/src/framework/containers/fields.h b/src/framework/containers/fields.h index ee9d656d6..4acabf7b4 100644 --- a/src/framework/containers/fields.h +++ b/src/framework/containers/fields.h @@ -21,6 +21,10 @@ #include "arch/kokkos_aliases.h" +#if defined(OUTPUT_ENABLED) + #include +#endif + #include namespace ntt { @@ -161,6 +165,18 @@ namespace ntt { (em_footprint + bckp_footprint + cur_footprint + buff_footprint + aux_footprint + em0_footprint + cur0_footprint); } + +/* helpers ---------------------------------------------------------------- */ +#if defined(OUTPUT_ENABLED) + void CheckpointDeclare(adios2::IO&, + const std::vector&, + const std::vector&, + const std::vector&) const; + void CheckpointRead(adios2::IO&, + adios2::Engine&, + const adios2::Box&); + void CheckpointWrite(adios2::IO&, adios2::Engine&) const; +#endif }; } // namespace ntt diff --git a/src/framework/containers/fields_io.cpp b/src/framework/containers/fields_io.cpp new file mode 100644 index 000000000..7d01451bd --- /dev/null +++ b/src/framework/containers/fields_io.cpp @@ -0,0 +1,96 @@ +#include "enums.h" +#include "global.h" + +#include "utils/log.h" + +#include "framework/containers/fields.h" +#include "output/utils/readers.h" +#include "output/utils/writers.h" + +#include + +#if defined(MPI_ENABLED) + #include +#endif + +#include + +namespace ntt { + + template + void Fields::CheckpointDeclare( + adios2::IO& io, + const std::vector& local_shape, + const std::vector& global_shape, + const std::vector& local_offset) const { + logger::Checkpoint("Declaring fields checkpoint", HERE); + + auto gs6 = std::vector(global_shape.begin(), global_shape.end()); + auto lo6 = std::vector(local_offset.begin(), local_offset.end()); + auto ls6 = std::vector(local_shape.begin(), local_shape.end()); + gs6.push_back(6); + lo6.push_back(0); + ls6.push_back(6); + + io.DefineVariable("em", gs6, lo6, ls6); + if (S == ntt::SimEngine::GRPIC) { + io.DefineVariable("em0", gs6, lo6, ls6); + auto gs3 = std::vector(global_shape.begin(), global_shape.end()); + auto lo3 = std::vector(local_offset.begin(), local_offset.end()); + auto ls3 = std::vector(local_shape.begin(), local_shape.end()); + gs3.push_back(3); + lo3.push_back(0); + ls3.push_back(3); + io.DefineVariable("cur0", gs3, lo3, ls3); + } + } + + template + void Fields::CheckpointRead(adios2::IO& io, + adios2::Engine& reader, + const adios2::Box& range) { + logger::Checkpoint("Reading fields checkpoint", HERE); + + auto range6 = adios2::Box(range.first, range.second); + range6.first.push_back(0); + range6.second.push_back(6); + out::ReadNDField(io, reader, "em", em, range6); + if (S == ntt::SimEngine::GRPIC) { + out::ReadNDField(io, reader, "em0", em0, range6); + auto range3 = adios2::Box(range.first, range.second); + range3.first.push_back(0); + range3.second.push_back(3); + out::ReadNDField(io, reader, "cur0", cur0, range3); + } + } + + template + void Fields::CheckpointWrite(adios2::IO& io, adios2::Engine& writer) const { + logger::Checkpoint("Writing fields checkpoint", HERE); + + out::WriteNDField(io, writer, "em", em); + if (S == ntt::SimEngine::GRPIC) { + out::WriteNDField(io, writer, "em0", em0); + out::WriteNDField(io, writer, "cur0", cur0); + } + } + +#define FIELDS_CHECKPOINTS(D, S) \ + template void Fields::CheckpointDeclare(adios2::IO&, \ + const std::vector&, \ + const std::vector&, \ + const std::vector&) \ + const; \ + template void Fields::CheckpointRead(adios2::IO&, \ + adios2::Engine&, \ + const adios2::Box&); \ + template void Fields::CheckpointWrite(adios2::IO&, adios2::Engine&) const; + + FIELDS_CHECKPOINTS(Dim::_1D, SimEngine::SRPIC) + FIELDS_CHECKPOINTS(Dim::_2D, SimEngine::SRPIC) + FIELDS_CHECKPOINTS(Dim::_3D, SimEngine::SRPIC) + FIELDS_CHECKPOINTS(Dim::_2D, SimEngine::GRPIC) + FIELDS_CHECKPOINTS(Dim::_3D, SimEngine::GRPIC) +#undef FIELDS_CHECKPOINTS + +} // namespace ntt diff --git a/src/framework/containers/particles_io.cpp b/src/framework/containers/particles_io.cpp index 870a6e0ae..70b7f2458 100644 --- a/src/framework/containers/particles_io.cpp +++ b/src/framework/containers/particles_io.cpp @@ -22,6 +22,7 @@ namespace ntt { logger::Checkpoint( fmt::format("Declaring particle checkpoint for species #%d", index()), HERE); + io.DefineVariable(fmt::format("s%d_npart", index()), { adios2::UnknownDim }, { adios2::UnknownDim }, @@ -97,24 +98,20 @@ namespace ntt { "Particles already initialized before reading checkpoint", HERE); npart_t npart_offset = 0u; + npart_t npart_read; out::ReadVariable(io, reader, fmt::format("s%d_npart", index()), - m_npart, + npart_read, domains_offset); - - raise::ErrorIf( - npart() > maxnpart(), - fmt::format("npart %d > maxnpart %d after reading checkpoint", - npart(), - maxnpart()), - HERE); + set_npart(npart_read); #if defined(MPI_ENABLED) { + const auto npart_send = npart(); std::vector glob_nparts(domains_total); - MPI_Allgather(&m_npart, + MPI_Allgather(&npart_send, 1, mpi::get_type(), glob_nparts.data(), @@ -324,24 +321,28 @@ namespace ntt { writer, fmt::format("s%d_i1", index()), i1, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_dx1", index()), dx1, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_i1_prev", index()), i1_prev, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_dx1_prev", index()), dx1_prev, + npart(), npart_total, npart_offset); } @@ -351,24 +352,28 @@ namespace ntt { writer, fmt::format("s%d_i2", index()), i2, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_dx2", index()), dx2, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_i2_prev", index()), i2_prev, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_dx2_prev", index()), dx2_prev, + npart(), npart_total, npart_offset); } @@ -378,24 +383,28 @@ namespace ntt { writer, fmt::format("s%d_i3", index()), i3, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_dx3", index()), dx3, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_i3_prev", index()), i3_prev, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_dx3_prev", index()), dx3_prev, + npart(), npart_total, npart_offset); } @@ -405,6 +414,7 @@ namespace ntt { writer, fmt::format("s%d_phi", index()), phi, + npart(), npart_total, npart_offset); } @@ -413,30 +423,35 @@ namespace ntt { writer, fmt::format("s%d_ux1", index()), ux1, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_ux2", index()), ux2, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_ux3", index()), ux3, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_tag", index()), tag, + npart(), npart_total, npart_offset); out::Write1DArray(io, writer, fmt::format("s%d_weight", index()), weight, + npart(), npart_total, npart_offset); if (npld_r() > 0) { @@ -444,9 +459,10 @@ namespace ntt { writer, fmt::format("s%d_pld_r", index()), pld_r, + npld_r(), + npart(), npart_total, - npart_offset, - npld_r()); + npart_offset); } if (npld_i() > 0) { @@ -454,29 +470,31 @@ namespace ntt { writer, fmt::format("s%d_pld_i", index()), pld_i, + npld_i(), + npart(), npart_total, - npart_offset, - npld_i()); + npart_offset); } } #define PARTICLES_CHECKPOINTS(D, C) \ - template void Particles::CheckpointDeclare(adios2::ADIOS&) const; \ - template void Particles::CheckpointRead(adios2::ADIOS&, \ + template void Particles::CheckpointDeclare(adios2::IO&) const; \ + template void Particles::CheckpointRead(adios2::IO&, \ adios2::Engine&, \ std::size_t, \ std::size_t); \ template void Particles::CheckpointWrite(adios2::IO&, \ adios2::Engine&, \ std::size_t, \ - std::size_t) const; \ - PARTICLES_CHECKPOINTS(Dim::_1D, Coord::Cart) \ - PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Cart) \ - PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Cart) \ - PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Sph) \ - PARTICLES_CHECKPOINTS(Dim::_2D, Coord::QSph) \ - PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Sph) \ - PARTICLES_CHECKPOINTS(Dim::_3D, Coord::QSph) + std::size_t) const; + + PARTICLES_CHECKPOINTS(Dim::_1D, Coord::Cart) + PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Cart) + PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Cart) + PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Sph) + PARTICLES_CHECKPOINTS(Dim::_2D, Coord::Qsph) + PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Sph) + PARTICLES_CHECKPOINTS(Dim::_3D, Coord::Qsph) #undef PARTICLES_CHECKPOINTS } // namespace ntt diff --git a/src/framework/domain/checkpoint.cpp b/src/framework/domain/checkpoint.cpp index 8fcefa989..c0a0fac6a 100644 --- a/src/framework/domain/checkpoint.cpp +++ b/src/framework/domain/checkpoint.cpp @@ -1,3 +1,5 @@ +#include "output/checkpoint.h" + #include "enums.h" #include "global.h" @@ -12,8 +14,6 @@ #include "metrics/qspherical.h" #include "metrics/spherical.h" -#include "checkpoint/reader.h" -#include "checkpoint/writer.h" #include "framework/domain/metadomain.h" #include "framework/parameters.h" @@ -59,11 +59,11 @@ namespace ntt { params.template get("checkpoint.keep"), params.template get("checkpoint.walltime")); if (g_checkpoint_writer.enabled()) { - g_checkpoint_writer.defineFieldVariables(S, - glob_shape_with_ghosts, - off_ncells_with_ghosts, - loc_shape_with_ghosts); - for (auto& species : local_domain->species) { + local_domain->fields.CheckpointDeclare(g_checkpoint_writer.io(), + loc_shape_with_ghosts, + glob_shape_with_ghosts, + off_ncells_with_ghosts); + for (const auto& species : local_domain->species) { species.CheckpointDeclare(g_checkpoint_writer.io()); } } @@ -90,16 +90,17 @@ namespace ntt { logger::Checkpoint("Writing checkpoint", HERE); g_checkpoint_writer.beginSaving(current_step, current_time); { - g_checkpoint_writer.saveAttrs(params, current_time); - g_checkpoint_writer.saveField("em", local_domain->fields.em); - if constexpr (S == SimEngine::GRPIC) { - g_checkpoint_writer.saveField("em0", local_domain->fields.em0); - g_checkpoint_writer.saveField("cur0", local_domain->fields.cur0); + if (g_checkpoint_writer.written().empty()) { + raise::Fatal("No checkpoint file to save metadata", HERE); } - std::size_t dom_tot = 1, dom_offset = 0; -#if defined(MPI_ENABLED) - dom_tot = g_mpi_size; - dom_offset = g_mpi_rank; + params.saveTOML(g_checkpoint_writer.written().back().second, current_time); + + local_domain->fields.CheckpointWrite(g_checkpoint_writer.io(), + g_checkpoint_writer.writer()); +#if !defined(MPI_ENABLED) + const std::size_t dom_tot = 1, dom_offset = 0; +#else + const std::size_t dom_tot = g_mpi_size, dom_offset = g_mpi_rank; #endif // MPI_ENABLED for (const auto& species : local_domain->species) { @@ -124,6 +125,7 @@ namespace ntt { fmt::format("step-%08lu.bp", params.template get( "checkpoint.start_step")); + logger::Checkpoint(fmt::format("Reading checkpoint from %s", fname.c_str()), HERE); @@ -137,45 +139,25 @@ namespace ntt { reader.BeginStep(); for (const auto local_domain_idx : l_subdomain_indices()) { - auto& domain = g_subdomains[local_domain_idx]; + auto local_domain = subdomain_ptr(local_domain_idx); + adios2::Box range; for (auto d { 0u }; d < M::Dim; ++d) { - range.first.push_back(domain.offset_ncells()[d] + - 2 * N_GHOSTS * domain.offset_ndomains()[d]); - range.second.push_back(domain.mesh.n_all()[d]); - } - range.first.push_back(0); - range.second.push_back(6); - checkpoint::ReadFields(io, reader, "em", range, domain.fields.em); - if constexpr (S == ntt::SimEngine::GRPIC) { - checkpoint::ReadFields(io, - reader, - "em0", - range, - domain.fields.em0); - adios2::Box range3; - for (auto d { 0u }; d < M::Dim; ++d) { - range3.first.push_back(domain.offset_ncells()[d] + - 2 * N_GHOSTS * domain.offset_ndomains()[d]); - range3.second.push_back(domain.mesh.n_all()[d]); - } - range3.first.push_back(0); - range3.second.push_back(3); - checkpoint::ReadFields(io, - reader, - "cur0", - range3, - domain.fields.cur0); + range.first.push_back(local_domain->offset_ncells()[d] + + 2 * N_GHOSTS * local_domain->offset_ndomains()[d]); + range.second.push_back(local_domain->mesh.n_all()[d]); } + local_domain->fields.CheckpointRead(io, reader, range); - for (auto& species : domain.species) { - species.CheckpointRead(io, reader, local_domain_idx, ndomains()); + for (auto& species : local_domain->species) { + species.CheckpointRead(io, reader, ndomains(), local_domain_idx); } } // local subdomain loop reader.EndStep(); reader.Close(); + logger::Checkpoint( fmt::format("Checkpoint reading done from %s", fname.c_str()), HERE); diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index 7ddacffb3..a456879b4 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -31,7 +31,7 @@ #endif // MPI_ENABLED #if defined(OUTPUT_ENABLED) - #include "checkpoint/writer.h" + #include "output/checkpoint.h" #include "output/writer.h" #include diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 27906b6e4..9ed1e1e18 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -23,6 +23,8 @@ #include #endif +#include +#include #include #include #include @@ -31,10 +33,10 @@ namespace ntt { template - auto get_dx0_V0(const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) - -> std::pair { + auto get_dx0_V0( + const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) -> std::pair { const auto metric = M(resolution, extent, params); const auto dx0 = metric.dxMin(); coord_t x_corner { ZERO }; @@ -1031,4 +1033,16 @@ namespace ntt { "Have not defined all the necessary variables", HERE); } + + void SimulationParams::saveTOML(const std::string& path, simtime_t time) const { + CallOnce([&]() { + std::ofstream metadata; + metadata.open(path); + metadata << "[metadata]\n" + << " time = " << time << "\n\n" + << data() << std::endl; + metadata.close(); + }); + } + } // namespace ntt diff --git a/src/framework/parameters.h b/src/framework/parameters.h index 723e860db..0af4fa405 100644 --- a/src/framework/parameters.h +++ b/src/framework/parameters.h @@ -20,6 +20,8 @@ #include "utils/param_container.h" #include "utils/toml.h" +#include + namespace ntt { struct SimulationParams : public prm::Parameters { @@ -52,6 +54,8 @@ namespace ntt { raw_data = data; } + void saveTOML(const std::string&, simtime_t) const; + private: toml::value raw_data; }; diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt index ef20dae56..e15378c62 100644 --- a/src/output/CMakeLists.txt +++ b/src/output/CMakeLists.txt @@ -30,6 +30,7 @@ set(SOURCES ${SRC_DIR}/stats.cpp ${SRC_DIR}/fields.cpp ${SRC_DIR}/utils/interpret_prompt.cpp) if(${output}) list(APPEND SOURCES ${SRC_DIR}/writer.cpp) + list(APPEND SOURCES ${SRC_DIR}/checkpoint.cpp) list(APPEND SOURCES ${SRC_DIR}/utils/writers.cpp) list(APPEND SOURCES ${SRC_DIR}/utils/readers.cpp) endif() diff --git a/src/output/checkpoint.cpp b/src/output/checkpoint.cpp new file mode 100644 index 000000000..4d74ae23e --- /dev/null +++ b/src/output/checkpoint.cpp @@ -0,0 +1,104 @@ +#include "output/checkpoint.h" + +#include "global.h" + +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/log.h" + +#include +#include + +#include +#include + +namespace checkpoint { + + void Writer::init(adios2::ADIOS* ptr_adios, + const path_t& checkpoint_root, + timestep_t interval, + simtime_t interval_time, + int keep, + const std::string& walltime) { + m_keep = keep; + m_checkpoint_root = checkpoint_root; + m_enabled = keep != 0; + if (not m_enabled) { + return; + } + m_tracker.init("checkpoint", interval, interval_time, walltime); + p_adios = ptr_adios; + raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE); + + m_io = p_adios->DeclareIO("Entity::Checkpoint"); + m_io.SetEngine("BPFile"); + + m_io.DefineVariable("Step"); + m_io.DefineVariable("Time"); + m_io.DefineAttribute("NGhosts", ntt::N_GHOSTS); + + CallOnce( + [](auto&& checkpoint_root) { + if (!std::filesystem::exists(checkpoint_root)) { + std::filesystem::create_directory(checkpoint_root); + } + }, + m_checkpoint_root); + } + + auto Writer::shouldSave(timestep_t step, simtime_t time) -> bool { + return m_enabled and m_tracker.shouldWrite(step, time); + } + + void Writer::beginSaving(timestep_t step, simtime_t time) { + raise::ErrorIf(!m_enabled, "Checkpoint is not enabled", HERE); + raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE); + if (m_writing_mode) { + raise::Fatal("Already writing", HERE); + } + m_writing_mode = true; + try { + const auto filename = m_checkpoint_root / fmt::format("step-%08lu.bp", step); + const auto metafilename = m_checkpoint_root / + fmt::format("meta-%08lu.toml", step); + m_writer = m_io.Open(filename, adios2::Mode::Write); + m_written.push_back({ filename, metafilename }); + logger::Checkpoint(fmt::format("Writing checkpoint to %s and %s", + filename.c_str(), + metafilename.c_str()), + HERE); + } catch (std::exception& e) { + raise::Fatal(e.what(), HERE); + } + + m_writer.BeginStep(); + m_writer.Put(m_io.InquireVariable("Step"), &step); + m_writer.Put(m_io.InquireVariable("Time"), &time); + } + + void Writer::endSaving() { + raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE); + if (!m_writing_mode) { + raise::Fatal("Not writing", HERE); + } + m_writing_mode = false; + m_writer.EndStep(); + m_writer.Close(); + + // optionally remove the oldest checkpoint + CallOnce([&]() { + if (m_keep > 0 and m_written.size() > (std::size_t)m_keep) { + const auto oldest = m_written.front(); + if (std::filesystem::exists(oldest.first) and + std::filesystem::exists(oldest.second)) { + std::filesystem::remove_all(oldest.first); + std::filesystem::remove(oldest.second); + m_written.erase(m_written.begin()); + } else { + raise::Warning("Checkpoint file does not exist for some reason", HERE); + } + } + }); + } + +} // namespace checkpoint diff --git a/src/checkpoint/writer.h b/src/output/checkpoint.h similarity index 65% rename from src/checkpoint/writer.h rename to src/output/checkpoint.h index 2750e2226..495b97b73 100644 --- a/src/checkpoint/writer.h +++ b/src/output/checkpoint.h @@ -1,24 +1,21 @@ /** - * @file checkpoint/writer.h - * @brief Class that dumps checkpoints + * @file output/checkpoint.h + * @brief Class that handles checkpoint writing * @implements * - checkpoint::Writer * @cpp: - * - writer.cpp + * - checkpoint.cpp * @namespaces: * - checkpoint:: */ -#ifndef CHECKPOINT_WRITER_H -#define CHECKPOINT_WRITER_H +#ifndef OUTPUT_CHECKPOINT_H +#define OUTPUT_CHECKPOINT_H -#include "enums.h" #include "global.h" #include "utils/tools.h" -#include "framework/parameters.h" - #include #include @@ -60,16 +57,6 @@ namespace checkpoint { void beginSaving(timestep_t, simtime_t); void endSaving(); - void saveAttrs(const ntt::SimulationParams&, simtime_t); - - template - void saveField(const std::string&, const ndfield_t&); - - void defineFieldVariables(const ntt::SimEngine&, - const std::vector&, - const std::vector&, - const std::vector&); - [[nodiscard]] auto io() -> adios2::IO& { return m_io; @@ -80,6 +67,11 @@ namespace checkpoint { return m_writer; } + [[nodiscard]] + auto written() const -> const std::vector>& { + return m_written; + } + [[nodiscard]] auto enabled() const -> bool { return m_enabled; @@ -88,4 +80,4 @@ namespace checkpoint { } // namespace checkpoint -#endif // CHECKPOINT_WRITER_H +#endif // OUTPUT_CHECKPOINT_H From 99137f41fbaf78c1c0ccee32a090be87009c8c30 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 21 Oct 2025 16:08:46 -0400 Subject: [PATCH 04/15] comm for particles placed in struct + pld comm --- src/engines/engine_printer.cpp | 3 +- src/framework/containers/particles.h | 19 + src/framework/containers/particles_comm.cpp | 394 ++++++++++++++++++++ src/framework/domain/comm_mpi.hpp | 289 -------------- src/framework/domain/comm_nompi.hpp | 2 - src/framework/domain/communications.cpp | 91 +---- src/kernels/comm.hpp | 74 ++-- 7 files changed, 485 insertions(+), 387 deletions(-) create mode 100644 src/framework/containers/particles_comm.cpp diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp index 20f5e81ba..cbfde9304 100644 --- a/src/engines/engine_printer.cpp +++ b/src/engines/engine_printer.cpp @@ -400,7 +400,8 @@ namespace ntt { add_param(report, 6, "GCA", "%s", species.use_gca() ? "ON" : "OFF"); } add_param(report, 6, "Cooling", "%s", species.cooling().to_string()); - add_param(report, 6, "# of payloads", "%d", species.npld()); + add_param(report, 6, "# of real-value payloads", "%d", species.npld_r()); + add_param(report, 6, "# of integer-value payloads", "%d", species.npld_i()); } report.pop_back(); }, diff --git a/src/framework/containers/particles.h b/src/framework/containers/particles.h index c85759035..47758b0d9 100644 --- a/src/framework/containers/particles.h +++ b/src/framework/containers/particles.h @@ -15,6 +15,7 @@ #include "enums.h" #include "global.h" +#include "arch/directions.h" #include "arch/kokkos_aliases.h" #include "utils/error.h" #include "utils/formatting.h" @@ -253,6 +254,24 @@ namespace ntt { */ void SyncHostDevice(); +#if defined(MPI_ENABLED) + /** + * @brief Communicate particles across neighboring meshblocks + * @param dirs_to_comm The directions requiring communication + * @param shifts_in_x1 The coordinate shifts in x1 direction per each communicated particle + * @param shifts_in_x2 The coordinate shifts in x2 direction per each communicated particle + * @param shifts_in_x3 The coordinate shifts in x3 direction per each communicated particle + * @param send_ranks The map of ranks per each send direction + * @param recv_ranks The map of ranks per each recv direction + */ + void Communicate(const dir::dirs_t&, + const array_t&, + const array_t&, + const array_t&, + const dir::map_t&, + const dir::map_t&); +#endif + #if defined(OUTPUT_ENABLED) void CheckpointDeclare(adios2::IO&) const; void CheckpointRead(adios2::IO&, adios2::Engine&, std::size_t, std::size_t); diff --git a/src/framework/containers/particles_comm.cpp b/src/framework/containers/particles_comm.cpp new file mode 100644 index 000000000..4d6d67118 --- /dev/null +++ b/src/framework/containers/particles_comm.cpp @@ -0,0 +1,394 @@ +#include "enums.h" +#include "global.h" + +#include "arch/directions.h" +#include "arch/kokkos_aliases.h" +#include "arch/mpi_aliases.h" +#include "arch/mpi_tags.h" +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/log.h" + +#include "framework/containers/particles.h" + +#include "kernels/comm.hpp" + +#include + +#include +#include + +namespace ntt { + + namespace prtls { + template + void send_recv(array_t& send_arr, + array_t& recv_arr, + int send_rank, + int recv_rank, + npart_t nsend, + npart_t nrecv, + npart_t offset) { +#if !defined(DEVICE_ENABLED) || defined(GPU_AWARE_MPI) + MPI_Sendrecv(send_arr.data(), + nsend, + mpi::get_type(), + send_rank, + 0, + recv_arr.data() + offset, + nrecv, + mpi::get_type(), + recv_rank, + 0, + MPI_COMM_WORLD, + MPI_STATUS_IGNORE); +#else + const auto slice = std::make_pair(offset, offset + nrecv); + + auto send_arr_h = Kokkos::create_mirror_view(send_arr); + auto recv_arr_h = Kokkos::create_mirror_view( + Kokkos::subview(recv_arr, slice)); + Kokkos::deep_copy(send_arr_h, send_arr); + MPI_Sendrecv(send_arr_h.data(), + nsend, + mpi::get_type(), + send_rank, + 0, + recv_arr_h.data(), + nrecv, + mpi::get_type(), + recv_rank, + 0, + MPI_COMM_WORLD, + MPI_STATUS_IGNORE); + Kokkos::deep_copy(Kokkos::subview(recv_arr, slice), recv_arr_h); +#endif + } + + void send_recv_count(int send_rank, + int recv_rank, + npart_t send_count, + npart_t& recv_count) { + if ((send_rank >= 0) && (recv_rank >= 0)) { + MPI_Sendrecv(&send_count, + 1, + mpi::get_type(), + send_rank, + 0, + &recv_count, + 1, + mpi::get_type(), + recv_rank, + 0, + MPI_COMM_WORLD, + MPI_STATUS_IGNORE); + } else if (send_rank >= 0) { + MPI_Send(&send_count, 1, mpi::get_type(), send_rank, 0, MPI_COMM_WORLD); + } else if (recv_rank >= 0) { + MPI_Recv(&recv_count, + 1, + mpi::get_type(), + recv_rank, + 0, + MPI_COMM_WORLD, + MPI_STATUS_IGNORE); + } else { + raise::Error("ParticleSendRecvCount called with negative ranks", HERE); + } + } + + template + void send(array_t& send_arr, int send_rank, npart_t nsend) { +#if !defined(DEVICE_ENABLED) || defined(GPU_AWARE_MPI) + MPI_Send(send_arr.data(), nsend, mpi::get_type(), send_rank, 0, MPI_COMM_WORLD); +#else + auto send_arr_h = Kokkos::create_mirror_view(send_arr); + Kokkos::deep_copy(send_arr_h, send_arr); + MPI_Send(send_arr_h.data(), nsend, mpi::get_type(), send_rank, 0, MPI_COMM_WORLD); +#endif + } + + template + void recv(array_t& recv_arr, int recv_rank, npart_t nrecv, npart_t offset) { +#if !defined(DEVICE_ENABLED) || defined(GPU_AWARE_MPI) + MPI_Recv(recv_arr.data() + offset, + nrecv, + mpi::get_type(), + recv_rank, + 0, + MPI_COMM_WORLD, + MPI_STATUS_IGNORE); +#else + const auto slice = std::make_pair(offset, offset + nrecv); + + auto recv_arr_h = Kokkos::create_mirror_view( + Kokkos::subview(recv_arr, slice)); + MPI_Recv(recv_arr_h.data(), + nrecv, + mpi::get_type(), + recv_rank, + 0, + MPI_COMM_WORLD, + MPI_STATUS_IGNORE); + Kokkos::deep_copy(Kokkos::subview(recv_arr, slice), recv_arr_h); +#endif + } + + template + void communicate(array_t& send_arr, + array_t& recv_arr, + int send_rank, + int recv_rank, + npart_t nsend, + npart_t nrecv, + npart_t offset) { + if (send_rank >= 0 && recv_rank >= 0) { + raise::ErrorIf( + nrecv + offset > recv_arr.extent(0), + "recv_arr is not large enough to hold the received particles", + HERE); + send_recv(send_arr, recv_arr, send_rank, recv_rank, nsend, nrecv, offset); + } else if (send_rank >= 0) { + send(send_arr, send_rank, nsend); + } else if (recv_rank >= 0) { + raise::ErrorIf( + nrecv + offset > recv_arr.extent(0), + "recv_arr is not large enough to hold the received particles", + HERE); + recv(recv_arr, recv_rank, nrecv, offset); + } else { + raise::Error("CommunicateParticles called with negative ranks", HERE); + } + } + } // namespace prtls + + template + void Particles::Communicate(const dir::dirs_t& dirs_to_comm, + const array_t& shifts_in_x1, + const array_t& shifts_in_x2, + const array_t& shifts_in_x3, + const dir::map_t& send_ranks, + const dir::map_t& recv_ranks) { + logger::Checkpoint(fmt::format("Communicating species #%d\n", index()), HERE); + + // at this point particles should already be tagged in the pusher + auto [npptag_vec, tag_offsets] = NpartsPerTagAndOffsets(); + const auto npart_dead = npptag_vec[ParticleTag::dead]; + const auto npart_alive = npptag_vec[ParticleTag::alive]; + + // # of particles to receive per each tag (direction) + std::vector npptag_recv_vec(ntags() - 2, 0); + + // total # of received particles from all directions + npart_t npart_recv_tot = 0u; + + // loop dir + for (const auto& direction : dirs_to_comm) { + // tags corresponding to the direction (both send & recv) + const auto tag_recv = mpi::PrtlSendTag::dir2tag(-direction); + const auto tag_send = mpi::PrtlSendTag::dir2tag(direction); + + // get ranks of send/recv meshblocks + const auto send_rank = send_ranks.at(direction); + const auto recv_rank = recv_ranks.at(direction); + + // record the # of particles to-be-sent + const auto nsend = npptag_vec[tag_send]; + + // request the # of particles to-be-received ... + // ... and send the # of particles to-be-sent + npart_t nrecv = 0; + prtls::send_recv_count(send_rank, recv_rank, nsend, nrecv); + npart_recv_tot += nrecv; + npptag_recv_vec[tag_recv - 2] = nrecv; + + raise::ErrorIf((npart() + npart_recv_tot) >= maxnpart(), + "Too many particles to receive (cannot fit into maxptl)", + HERE); + } + + array_t outgoing_indices { "outgoing_indices", npart() - npart_alive }; + // clang-format off + Kokkos::parallel_for( + "PrepareOutgoingPrtls", + rangeActiveParticles(), + kernel::comm::PrepareOutgoingPrtls_kernel( + shifts_in_x1, shifts_in_x2, shifts_in_x3, + outgoing_indices, + npart(), npart_alive, npart_dead, ntags(), + i1, i1_prev, + i2, i2_prev, + i3, i3_prev, + tag, tag_offsets) + ); + // clang-format on + + // number of arrays of each type to send/recv + const unsigned short NREALS = 4 + static_cast( + D == Dim::_2D and C != Coord::Cart); + const unsigned short NINTS = 2 * static_cast(D); + const unsigned short NPRTLDX = 2 * static_cast(D); + const unsigned short NPLDS_R = npld_r(); + const unsigned short NPLDS_I = npld_i(); + + // buffers to store recv data + const auto npart_recv = std::accumulate(npptag_recv_vec.begin(), + npptag_recv_vec.end(), + static_cast(0)); + array_t recv_buff_int { "recv_buff_int", npart_recv * NINTS }; + array_t recv_buff_real { "recv_buff_real", npart_recv * NREALS }; + array_t recv_buff_prtldx { "recv_buff_prtldx", npart_recv * NPRTLDX }; + array_t recv_buff_pld_r; + array_t recv_buff_pld_i; + + if (NPLDS_R > 0) { + recv_buff_pld_r = array_t { "recv_buff_pld_r", npart_recv * NPLDS_R }; + } + if (NPLDS_I > 0) { + recv_buff_pld_i = array_t { "recv_buff_pld_i", + npart_recv * NPLDS_I }; + } + + auto iteration = 0; + auto current_received = 0; + + for (const auto& direction : dirs_to_comm) { + const auto send_rank = send_ranks.at(direction); + const auto recv_rank = recv_ranks.at(direction); + const auto tag_send = mpi::PrtlSendTag::dir2tag(direction); + const auto tag_recv = mpi::PrtlSendTag::dir2tag(-direction); + const auto npart_send_in = npptag_vec[tag_send]; + const auto npart_recv_in = npptag_recv_vec[tag_recv - 2]; + if (send_rank < 0 and recv_rank < 0) { + continue; + } + array_t send_buff_int { "send_buff_int", npart_send_in * NINTS }; + array_t send_buff_real { "send_buff_real", npart_send_in * NREALS }; + array_t send_buff_prtldx { "send_buff_prtldx", + npart_send_in * NPRTLDX }; + array_t send_buff_pld_r; + array_t send_buff_pld_i; + if (NPLDS_R > 0) { + send_buff_pld_r = array_t { "send_buff_pld_r", + npart_send_in * NPLDS_R }; + } + if (NPLDS_I > 0) { + send_buff_pld_i = array_t { "send_buff_pld_i", + npart_send_in * NPLDS_I }; + } + + auto tag_offsets_h = Kokkos::create_mirror_view(tag_offsets); + Kokkos::deep_copy(tag_offsets_h, tag_offsets); + + npart_t idx_offset = npart_dead; + if (tag_send > 2) { + idx_offset += tag_offsets_h(tag_send - 3); + } + // clang-format off + Kokkos::parallel_for( + "PopulatePrtlSendBuffer", + npart_send_in, + kernel::comm::PopulatePrtlSendBuffer_kernel( + send_buff_int, send_buff_real, send_buff_prtldx, send_buff_pld_r, send_buff_pld_i, + NINTS, NREALS, NPRTLDX, NPLDS_R, NPLDS_I, idx_offset, + i1, i1_prev, dx1, dx1_prev, + i2, i2_prev, dx2, dx2_prev, + i3, i3_prev, dx3, dx3_prev, + ux1, ux2, ux3, + weight, phi, pld_r, pld_i, tag, + outgoing_indices) + ); + // clang-format on + + const auto recv_offset_int = current_received * NINTS; + const auto recv_offset_real = current_received * NREALS; + const auto recv_offset_prtldx = current_received * NPRTLDX; + const auto recv_offset_pld_r = current_received * NPLDS_R; + const auto recv_offset_pld_i = current_received * NPLDS_I; + + prtls::communicate(send_buff_int, + recv_buff_int, + send_rank, + recv_rank, + npart_send_in * NINTS, + npart_recv_in * NINTS, + recv_offset_int); + prtls::communicate(send_buff_real, + recv_buff_real, + send_rank, + recv_rank, + npart_send_in * NREALS, + npart_recv_in * NREALS, + recv_offset_real); + prtls::communicate(send_buff_prtldx, + recv_buff_prtldx, + send_rank, + recv_rank, + npart_send_in * NPRTLDX, + npart_recv_in * NPRTLDX, + recv_offset_prtldx); + if (NPLDS_R > 0) { + prtls::communicate(send_buff_pld_r, + recv_buff_pld_r, + send_rank, + recv_rank, + npart_send_in * NPLDS_R, + npart_recv_in * NPLDS_R, + recv_offset_pld_r); + } + if (NPLDS_I > 0) { + prtls::communicate(send_buff_pld_i, + recv_buff_pld_i, + send_rank, + recv_rank, + npart_send_in * NPLDS_I, + npart_recv_in * NPLDS_I, + recv_offset_pld_i); + } + current_received += npart_recv_in; + iteration++; + + } // end direction loop + + // clang-format off + Kokkos::parallel_for( + "PopulateFromRecvBuffer", + npart_recv, + kernel::comm::ExtractReceivedPrtls_kernel( + recv_buff_int, recv_buff_real, recv_buff_prtldx, recv_buff_pld_r, recv_buff_pld_i, + NINTS, NREALS, NPRTLDX, NPLDS_R, NPLDS_I, + npart(), + i1, i1_prev, dx1, dx1_prev, + i2, i2_prev, dx2, dx2_prev, + i3, i3_prev, dx3, dx3_prev, + ux1, ux2, ux3, + weight, phi, pld_r, pld_i, tag, + outgoing_indices) + ); + // clang-format on + + const auto npart_holes = outgoing_indices.extent(0); + if (npart_recv > npart_holes) { + set_npart(npart() + npart_recv - npart_holes); + } + set_unsorted(); + } + +#define PARTICLES_COMM(D, C) \ + template void Particles::Communicate(const dir::dirs_t&, \ + const array_t&, \ + const array_t&, \ + const array_t&, \ + const dir::map_t&, \ + const dir::map_t&); + + PARTICLES_COMM(Dim::_1D, Coord::Cart) + PARTICLES_COMM(Dim::_2D, Coord::Cart) + PARTICLES_COMM(Dim::_3D, Coord::Cart) + PARTICLES_COMM(Dim::_2D, Coord::Sph) + PARTICLES_COMM(Dim::_2D, Coord::Qsph) + PARTICLES_COMM(Dim::_3D, Coord::Sph) + PARTICLES_COMM(Dim::_3D, Coord::Qsph) +#undef PARTICLES_COMM + +} // namespace ntt diff --git a/src/framework/domain/comm_mpi.hpp b/src/framework/domain/comm_mpi.hpp index e0d0cb4b2..266c844c5 100644 --- a/src/framework/domain/comm_mpi.hpp +++ b/src/framework/domain/comm_mpi.hpp @@ -11,23 +11,15 @@ #ifndef FRAMEWORK_DOMAIN_COMM_MPI_HPP #define FRAMEWORK_DOMAIN_COMM_MPI_HPP -#include "enums.h" #include "global.h" -#include "arch/directions.h" #include "arch/kokkos_aliases.h" #include "arch/mpi_aliases.h" -#include "arch/mpi_tags.h" #include "utils/error.h" -#include "framework/containers/particles.h" - -#include "kernels/comm.hpp" - #include #include -#include #include namespace comm { @@ -131,116 +123,6 @@ namespace comm { } // namespace flds - namespace prtls { - template - void send_recv(array_t& send_arr, - array_t& recv_arr, - int send_rank, - int recv_rank, - npart_t nsend, - npart_t nrecv, - npart_t offset) { -#if !defined(DEVICE_ENABLED) || defined(GPU_AWARE_MPI) - MPI_Sendrecv(send_arr.data(), - nsend, - mpi::get_type(), - send_rank, - 0, - recv_arr.data() + offset, - nrecv, - mpi::get_type(), - recv_rank, - 0, - MPI_COMM_WORLD, - MPI_STATUS_IGNORE); -#else - const auto slice = std::make_pair(offset, offset + nrecv); - - auto send_arr_h = Kokkos::create_mirror_view(send_arr); - auto recv_arr_h = Kokkos::create_mirror_view( - Kokkos::subview(recv_arr, slice)); - Kokkos::deep_copy(send_arr_h, send_arr); - MPI_Sendrecv(send_arr_h.data(), - nsend, - mpi::get_type(), - send_rank, - 0, - recv_arr_h.data(), - nrecv, - mpi::get_type(), - recv_rank, - 0, - MPI_COMM_WORLD, - MPI_STATUS_IGNORE); - Kokkos::deep_copy(Kokkos::subview(recv_arr, slice), recv_arr_h); -#endif - } - - template - void send(array_t& send_arr, int send_rank, npart_t nsend) { -#if !defined(DEVICE_ENABLED) || defined(GPU_AWARE_MPI) - MPI_Send(send_arr.data(), nsend, mpi::get_type(), send_rank, 0, MPI_COMM_WORLD); -#else - auto send_arr_h = Kokkos::create_mirror_view(send_arr); - Kokkos::deep_copy(send_arr_h, send_arr); - MPI_Send(send_arr_h.data(), nsend, mpi::get_type(), send_rank, 0, MPI_COMM_WORLD); -#endif - } - - template - void recv(array_t& recv_arr, int recv_rank, npart_t nrecv, npart_t offset) { -#if !defined(DEVICE_ENABLED) || defined(GPU_AWARE_MPI) - MPI_Recv(recv_arr.data() + offset, - nrecv, - mpi::get_type(), - recv_rank, - 0, - MPI_COMM_WORLD, - MPI_STATUS_IGNORE); -#else - const auto slice = std::make_pair(offset, offset + nrecv); - - auto recv_arr_h = Kokkos::create_mirror_view( - Kokkos::subview(recv_arr, slice)); - MPI_Recv(recv_arr_h.data(), - nrecv, - mpi::get_type(), - recv_rank, - 0, - MPI_COMM_WORLD, - MPI_STATUS_IGNORE); - Kokkos::deep_copy(Kokkos::subview(recv_arr, slice), recv_arr_h); -#endif - } - - template - void communicate(array_t& send_arr, - array_t& recv_arr, - int send_rank, - int recv_rank, - npart_t nsend, - npart_t nrecv, - npart_t offset) { - if (send_rank >= 0 && recv_rank >= 0) { - raise::ErrorIf( - nrecv + offset > recv_arr.extent(0), - "recv_arr is not large enough to hold the received particles", - HERE); - send_recv(send_arr, recv_arr, send_rank, recv_rank, nsend, nrecv, offset); - } else if (send_rank >= 0) { - send(send_arr, send_rank, nsend); - } else if (recv_rank >= 0) { - raise::ErrorIf( - nrecv + offset > recv_arr.extent(0), - "recv_arr is not large enough to hold the received particles", - HERE); - recv(recv_arr, recv_rank, nrecv, offset); - } else { - raise::Error("CommunicateParticles called with negative ranks", HERE); - } - } - } // namespace prtls - template inline void CommunicateField(unsigned int idx, ndfield_t& fld, @@ -468,177 +350,6 @@ namespace comm { } } - void ParticleSendRecvCount(int send_rank, - int recv_rank, - npart_t send_count, - npart_t& recv_count) { - if ((send_rank >= 0) && (recv_rank >= 0)) { - MPI_Sendrecv(&send_count, - 1, - mpi::get_type(), - send_rank, - 0, - &recv_count, - 1, - mpi::get_type(), - recv_rank, - 0, - MPI_COMM_WORLD, - MPI_STATUS_IGNORE); - } else if (send_rank >= 0) { - MPI_Send(&send_count, 1, mpi::get_type(), send_rank, 0, MPI_COMM_WORLD); - } else if (recv_rank >= 0) { - MPI_Recv(&recv_count, - 1, - mpi::get_type(), - recv_rank, - 0, - MPI_COMM_WORLD, - MPI_STATUS_IGNORE); - } else { - raise::Error("ParticleSendRecvCount called with negative ranks", HERE); - } - } - - template - void CommunicateParticles(Particles& species, - const array_t& outgoing_indices, - const array_t& tag_offsets, - const std::vector& npptag_vec, - const std::vector& npptag_recv_vec, - const std::vector& send_ranks, - const std::vector& recv_ranks, - const dir::dirs_t& dirs_to_comm) { - // number of arrays of each type to send/recv - const unsigned short NREALS = 4 + static_cast( - D == Dim::_2D and C != Coord::Cart); - const unsigned short NINTS = 2 * static_cast(D); - const unsigned short NPRTLDX = 2 * static_cast(D); - const unsigned short NPLDS = species.npld(); - - // buffers to store recv data - const auto npart_dead = npptag_vec[ParticleTag::dead]; - const auto npart_recv = std::accumulate(npptag_recv_vec.begin(), - npptag_recv_vec.end(), - static_cast(0)); - array_t recv_buff_int { "recv_buff_int", npart_recv * NINTS }; - array_t recv_buff_real { "recv_buff_real", npart_recv * NREALS }; - array_t recv_buff_prtldx { "recv_buff_prtldx", npart_recv * NPRTLDX }; - array_t recv_buff_pld; - - if (NPLDS > 0) { - recv_buff_pld = array_t { "recv_buff_pld", npart_recv * NPLDS }; - } - - auto iteration = 0; - auto current_received = 0; - - for (const auto& direction : dirs_to_comm) { - const auto send_rank = send_ranks[iteration]; - const auto recv_rank = recv_ranks[iteration]; - const auto tag_send = mpi::PrtlSendTag::dir2tag(direction); - const auto tag_recv = mpi::PrtlSendTag::dir2tag(-direction); - const auto npart_send_in = npptag_vec[tag_send]; - const auto npart_recv_in = npptag_recv_vec[tag_recv - 2]; - if (send_rank < 0 and recv_rank < 0) { - continue; - } - array_t send_buff_int { "send_buff_int", npart_send_in * NINTS }; - array_t send_buff_real { "send_buff_real", npart_send_in * NREALS }; - array_t send_buff_prtldx { "send_buff_prtldx", - npart_send_in * NPRTLDX }; - array_t send_buff_pld; - if (NPLDS > 0) { - send_buff_pld = array_t { "send_buff_pld", npart_send_in * NPLDS }; - } - - auto tag_offsets_h = Kokkos::create_mirror_view(tag_offsets); - Kokkos::deep_copy(tag_offsets_h, tag_offsets); - - npart_t idx_offset = npart_dead; - if (tag_send > 2) { - idx_offset += tag_offsets_h(tag_send - 3); - } - // clang-format off - Kokkos::parallel_for( - "PopulatePrtlSendBuffer", - npart_send_in, - kernel::comm::PopulatePrtlSendBuffer_kernel( - send_buff_int, send_buff_real, send_buff_prtldx, send_buff_pld, - NINTS, NREALS, NPRTLDX, NPLDS, idx_offset, - species.i1, species.i1_prev, species.dx1, species.dx1_prev, - species.i2, species.i2_prev, species.dx2, species.dx2_prev, - species.i3, species.i3_prev, species.dx3, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.weight, species.phi, species.pld, species.tag, - outgoing_indices) - ); - // clang-format on - - const auto recv_offset_int = current_received * NINTS; - const auto recv_offset_real = current_received * NREALS; - const auto recv_offset_prtldx = current_received * NPRTLDX; - const auto recv_offset_pld = current_received * NPLDS; - - prtls::communicate(send_buff_int, - recv_buff_int, - send_rank, - recv_rank, - npart_send_in * NINTS, - npart_recv_in * NINTS, - recv_offset_int); - prtls::communicate(send_buff_real, - recv_buff_real, - send_rank, - recv_rank, - npart_send_in * NREALS, - npart_recv_in * NREALS, - recv_offset_real); - prtls::communicate(send_buff_prtldx, - recv_buff_prtldx, - send_rank, - recv_rank, - npart_send_in * NPRTLDX, - npart_recv_in * NPRTLDX, - recv_offset_prtldx); - if (NPLDS > 0) { - prtls::communicate(send_buff_pld, - recv_buff_pld, - send_rank, - recv_rank, - npart_send_in * NPLDS, - npart_recv_in * NPLDS, - recv_offset_pld); - } - current_received += npart_recv_in; - iteration++; - - } // end direction loop - - // clang-format off - Kokkos::parallel_for( - "PopulateFromRecvBuffer", - npart_recv, - kernel::comm::ExtractReceivedPrtls_kernel( - recv_buff_int, recv_buff_real, recv_buff_prtldx, recv_buff_pld, - NINTS, NREALS, NPRTLDX, NPLDS, - species.npart(), - species.i1, species.i1_prev, species.dx1, species.dx1_prev, - species.i2, species.i2_prev, species.dx2, species.dx2_prev, - species.i3, species.i3_prev, species.dx3, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.weight, species.phi, species.pld, species.tag, - outgoing_indices) - ); - // clang-format on - - const auto npart = species.npart(); - const auto npart_holes = outgoing_indices.extent(0); - if (npart_recv > npart_holes) { - species.set_npart(npart + npart_recv - npart_holes); - } - } - } // namespace comm #endif // FRAMEWORK_DOMAIN_COMM_MPI_HPP diff --git a/src/framework/domain/comm_nompi.hpp b/src/framework/domain/comm_nompi.hpp index b477ac176..d0cc36cbd 100644 --- a/src/framework/domain/comm_nompi.hpp +++ b/src/framework/domain/comm_nompi.hpp @@ -16,8 +16,6 @@ #include "arch/kokkos_aliases.h" #include "utils/error.h" -#include "framework/domain/domain.h" - #include namespace comm { diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp index bf6eb3dd1..59ddbf583 100644 --- a/src/framework/domain/communications.cpp +++ b/src/framework/domain/communications.cpp @@ -5,7 +5,6 @@ #include "utils/error.h" #include "utils/formatting.h" #include "utils/log.h" -#include "utils/timer.h" #include "metrics/kerr_schild.h" #include "metrics/kerr_schild_0.h" @@ -20,7 +19,6 @@ #include "arch/mpi_tags.h" #include "framework/domain/comm_mpi.hpp" - #include "kernels/comm.hpp" #else #include "framework/domain/comm_nompi.hpp" #endif @@ -587,45 +585,30 @@ namespace ntt { for (auto& species : domain.species) { const auto ntags = species.ntags(); - // at this point particles should already be tagged in the pusher - auto [npptag_vec, tag_offsets] = species.NpartsPerTagAndOffsets(); - const auto npart_dead = npptag_vec[ParticleTag::dead]; - const auto npart_alive = npptag_vec[ParticleTag::alive]; - - const auto npart = species.npart(); - - // # of particles to receive per each tag (direction) - std::vector npptag_recv_vec(ntags - 2, 0); // coordinate shifts per each direction - array_t shifts_in_x1 { "shifts_in_x1", ntags - 2 }; - array_t shifts_in_x2 { "shifts_in_x2", ntags - 2 }; - array_t shifts_in_x3 { "shifts_in_x3", ntags - 2 }; - auto shifts_in_x1_h = Kokkos::create_mirror_view(shifts_in_x1); - auto shifts_in_x2_h = Kokkos::create_mirror_view(shifts_in_x2); - auto shifts_in_x3_h = Kokkos::create_mirror_view(shifts_in_x3); + array_t shifts_in_x1 { "shifts_in_x1", ntags - 2 }; + array_t shifts_in_x2 { "shifts_in_x2", ntags - 2 }; + array_t shifts_in_x3 { "shifts_in_x3", ntags - 2 }; + auto shifts_in_x1_h = Kokkos::create_mirror_view(shifts_in_x1); + auto shifts_in_x2_h = Kokkos::create_mirror_view(shifts_in_x2); + auto shifts_in_x3_h = Kokkos::create_mirror_view(shifts_in_x3); // all directions requiring communication dir::dirs_t dirs_to_comm; // ranks & indices of meshblock to send/recv from - std::vector send_ranks, send_inds; - std::vector recv_ranks, recv_inds; - - // total # of reaceived particles from all directions - npart_t npart_recv = 0u; + dir::map_t send_ranks; + dir::map_t recv_ranks; for (const auto& direction : dir::Directions::all) { // tags corresponding to the direction (both send & recv) - const auto tag_recv = mpi::PrtlSendTag::dir2tag(-direction); const auto tag_send = mpi::PrtlSendTag::dir2tag(direction); // get indices & ranks of send/recv meshblocks const auto [send_params, - recv_params] = GetSendRecvParams(this, domain, direction, true); - const auto [send_indrank, send_slice] = send_params; - const auto [recv_indrank, recv_slice] = recv_params; - const auto [send_ind, send_rank] = send_indrank; - const auto [recv_ind, recv_rank] = recv_indrank; + recv_params] = GetSendRecvRanks(this, domain, direction); + const auto [send_ind, send_rank] = send_params; + const auto [recv_ind, recv_rank] = recv_params; // skip if no communication is necessary const auto is_sending = (send_rank >= 0); @@ -634,24 +617,8 @@ namespace ntt { continue; } dirs_to_comm.push_back(direction); - send_ranks.push_back(send_rank); - recv_ranks.push_back(recv_rank); - send_inds.push_back(send_ind); - recv_inds.push_back(recv_ind); - - // record the # of particles to-be-sent - const auto nsend = npptag_vec[tag_send]; - - // request the # of particles to-be-received ... - // ... and send the # of particles to-be-sent - npart_t nrecv = 0; - comm::ParticleSendRecvCount(send_rank, recv_rank, nsend, nrecv); - npart_recv += nrecv; - npptag_recv_vec[tag_recv - 2] = nrecv; - - raise::ErrorIf((npart + npart_recv) >= species.maxnpart(), - "Too many particles to receive (cannot fit into maxptl)", - HERE); + send_ranks[direction] = send_rank; + recv_ranks[direction] = recv_rank; // if sending, record displacements to apply before // ... tag_send - 2: because we only shift tags > 2 (i.e. no dead/alive) @@ -689,31 +656,13 @@ namespace ntt { Kokkos::deep_copy(shifts_in_x2, shifts_in_x2_h); Kokkos::deep_copy(shifts_in_x3, shifts_in_x3_h); - array_t outgoing_indices { "outgoing_indices", npart - npart_alive }; - // clang-format off - Kokkos::parallel_for( - "PrepareOutgoingPrtls", - species.rangeActiveParticles(), - kernel::comm::PrepareOutgoingPrtls_kernel( - shifts_in_x1, shifts_in_x2, shifts_in_x3, - outgoing_indices, - npart, npart_alive, npart_dead, ntags, - species.i1, species.i1_prev, - species.i2, species.i2_prev, - species.i3, species.i3_prev, - species.tag, tag_offsets) - ); - // clang-format on - - comm::CommunicateParticles(species, - outgoing_indices, - tag_offsets, - npptag_vec, - npptag_recv_vec, - send_ranks, - recv_ranks, - dirs_to_comm); - species.set_unsorted(); + species.Communicate(dirs_to_comm, + shifts_in_x1, + shifts_in_x2, + shifts_in_x3, + send_ranks, + recv_ranks); + } // end species loop #else (void)domain; diff --git a/src/kernels/comm.hpp b/src/kernels/comm.hpp index 60251d8c6..9f2645cdc 100644 --- a/src/kernels/comm.hpp +++ b/src/kernels/comm.hpp @@ -110,15 +110,17 @@ namespace kernel::comm { array_t send_buff_int; array_t send_buff_real; array_t send_buff_prtldx; - array_t send_buff_pld; + array_t send_buff_pld_r; + array_t send_buff_pld_i; - const unsigned short NINTS, NREALS, NPRTLDX, NPLDS; + const unsigned short NINTS, NREALS, NPRTLDX, NPLDS_R, NPLDS_I; const npart_t idx_offset; const array_t i1, i1_prev, i2, i2_prev, i3, i3_prev; const array_t dx1, dx1_prev, dx2, dx2_prev, dx3, dx3_prev; const array_t ux1, ux2, ux3, weight, phi; - const array_t pld; + const array_t pld_r; + const array_t pld_i; array_t tag; const array_t outgoing_indices; @@ -126,11 +128,13 @@ namespace kernel::comm { PopulatePrtlSendBuffer_kernel(array_t& send_buff_int, array_t& send_buff_real, array_t& send_buff_prtldx, - array_t& send_buff_pld, + array_t& send_buff_pld_r, + array_t& send_buff_pld_i, unsigned short NINTS, unsigned short NREALS, unsigned short NPRTLDX, - unsigned short NPLDS, + unsigned short NPLDS_R, + unsigned short NPLDS_I, npart_t idx_offset, const array_t& i1, const array_t& i1_prev, @@ -149,17 +153,20 @@ namespace kernel::comm { const array_t& ux3, const array_t& weight, const array_t& phi, - const array_t& pld, + const array_t& pld_r, + const array_t& pld_i, array_t& tag, const array_t& outgoing_indices) : send_buff_int { send_buff_int } , send_buff_real { send_buff_real } , send_buff_prtldx { send_buff_prtldx } - , send_buff_pld { send_buff_pld } + , send_buff_pld_r { send_buff_pld_r } + , send_buff_pld_i { send_buff_pld_i } , NINTS { NINTS } , NREALS { NREALS } , NPRTLDX { NPRTLDX } - , NPLDS { NPLDS } + , NPLDS_R { NPLDS_R } + , NPLDS_I { NPLDS_I } , idx_offset { idx_offset } , i1 { i1 } , i1_prev { i1_prev } @@ -178,7 +185,8 @@ namespace kernel::comm { , ux3 { ux3 } , weight { weight } , phi { phi } - , pld { pld } + , pld_r { pld_r } + , pld_i { pld_i } , tag { tag } , outgoing_indices { outgoing_indices } {} @@ -209,9 +217,14 @@ namespace kernel::comm { if constexpr (D == Dim::_2D and C != Coord::Cart) { send_buff_real(NREALS * p + 4) = phi(idx); } - if (NPLDS > 0) { - for (auto l { 0u }; l < NPLDS; ++l) { - send_buff_pld(NPLDS * p + l) = pld(idx, l); + if (NPLDS_R > 0) { + for (auto l { 0u }; l < NPLDS_R; ++l) { + send_buff_pld_r(NPLDS_R * p + l) = pld_r(idx, l); + } + } + if (NPLDS_I > 0) { + for (auto l { 0u }; l < NPLDS_I; ++l) { + send_buff_pld_i(NPLDS_I * p + l) = pld_i(idx, l); } } tag(idx) = ParticleTag::dead; @@ -223,15 +236,17 @@ namespace kernel::comm { const array_t recv_buff_int; const array_t recv_buff_real; const array_t recv_buff_prtldx; - const array_t recv_buff_pld; + const array_t recv_buff_pld_r; + const array_t recv_buff_pld_i; - const unsigned short NINTS, NREALS, NPRTLDX, NPLDS; + const unsigned short NINTS, NREALS, NPRTLDX, NPLDS_R, NPLDS_I; const npart_t npart, npart_holes; array_t i1, i1_prev, i2, i2_prev, i3, i3_prev; array_t dx1, dx1_prev, dx2, dx2_prev, dx3, dx3_prev; array_t ux1, ux2, ux3, weight, phi; - array_t pld; + array_t pld_r; + array_t pld_i; array_t tag; const array_t outgoing_indices; @@ -239,11 +254,13 @@ namespace kernel::comm { ExtractReceivedPrtls_kernel(const array_t& recv_buff_int, const array_t& recv_buff_real, const array_t& recv_buff_prtldx, - const array_t& recv_buff_pld, + const array_t& recv_buff_pld_r, + const array_t& recv_buff_pld_i, unsigned short NINTS, unsigned short NREALS, unsigned short NPRTLDX, - unsigned short NPLDS, + unsigned short NPLDS_R, + unsigned short NPLDS_I, npart_t npart, array_t& i1, array_t& i1_prev, @@ -262,17 +279,20 @@ namespace kernel::comm { array_t& ux3, array_t& weight, array_t& phi, - array_t& pld, + array_t& pld_r, + array_t& pld_i, array_t& tag, const array_t& outgoing_indices) : recv_buff_int { recv_buff_int } , recv_buff_real { recv_buff_real } , recv_buff_prtldx { recv_buff_prtldx } - , recv_buff_pld { recv_buff_pld } + , recv_buff_pld_r { recv_buff_pld_r } + , recv_buff_pld_i { recv_buff_pld_i } , NINTS { NINTS } , NREALS { NREALS } , NPRTLDX { NPRTLDX } - , NPLDS { NPLDS } + , NPLDS_R { NPLDS_R } + , NPLDS_I { NPLDS_I } , npart { npart } , npart_holes { outgoing_indices.extent(0) } , i1 { i1 } @@ -292,7 +312,8 @@ namespace kernel::comm { , ux3 { ux3 } , weight { weight } , phi { phi } - , pld { pld } + , pld_r { pld_r } + , pld_i { pld_i } , tag { tag } , outgoing_indices { outgoing_indices } {} @@ -328,9 +349,14 @@ namespace kernel::comm { if constexpr (D == Dim::_2D and C != Coord::Cart) { phi(idx) = recv_buff_real(NREALS * p + 4); } - if (NPLDS > 0) { - for (auto l { 0u }; l < NPLDS; ++l) { - pld(idx, l) = recv_buff_pld(NPLDS * p + l); + if (NPLDS_R > 0) { + for (auto l { 0u }; l < NPLDS_R; ++l) { + pld_r(idx, l) = recv_buff_pld_r(NPLDS_R * p + l); + } + } + if (NPLDS_I > 0) { + for (auto l { 0u }; l < NPLDS_I; ++l) { + pld_i(idx, l) = recv_buff_pld_i(NPLDS_I * p + l); } } tag(idx) = ParticleTag::alive; From 53bff3c6cdb9b32983edf27f36dbd38638a4aacc Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 21 Oct 2025 16:08:58 -0400 Subject: [PATCH 05/15] generics for write ndfield --- src/output/utils/readers.cpp | 49 +++++++++++++++++++++++++++++++----- src/output/utils/readers.h | 8 ++++++ src/output/utils/writers.cpp | 46 +++++++++++++++++++++++++-------- src/output/utils/writers.h | 7 ++++++ 4 files changed, 94 insertions(+), 16 deletions(-) diff --git a/src/output/utils/readers.cpp b/src/output/utils/readers.cpp index d2e866af1..53746dcec 100644 --- a/src/output/utils/readers.cpp +++ b/src/output/utils/readers.cpp @@ -20,8 +20,10 @@ namespace out { std::size_t local_offset) { auto var = io.InquireVariable(quantity); if (var) { + T read_data; var.SetSelection(adios2::Box({ local_offset }, { 1 })); - reader.Get(var, &data, adios2::Mode::Sync); + reader.Get(var, &read_data, adios2::Mode::Sync); + data = read_data; } else { raise::Error(fmt::format("Variable: %s not found", quantity.c_str()), HERE); } @@ -70,6 +72,24 @@ namespace out { } } + template + void ReadNDField(adios2::IO& io, + adios2::Engine& reader, + const std::string& quantity, + ndfield_t& data, + const adios2::Box& range) { + auto var = io.InquireVariable(quantity); + if (var) { + var.SetSelection(range); + + auto data_h = Kokkos::create_mirror_view(data); + reader.Get(var, data_h.data(), adios2::Mode::Sync); + Kokkos::deep_copy(data, data_h); + } else { + raise::Error(fmt::format("Variable: %s not found", quantity.c_str()), HERE); + } + } + #define ARRAY_READERS(T) \ template void ReadVariable(adios2::IO&, \ adios2::Engine&, \ @@ -88,12 +108,29 @@ namespace out { array_t&, \ unsigned short, \ std::size_t, \ - std::size_t); \ - ARRAY_READERS(int) \ - ARRAY_READERS(unsigned int) \ - ARRAY_READERS(unsigned long int) \ - ARRAY_READERS(double) \ + std::size_t); + + ARRAY_READERS(short) + ARRAY_READERS(unsigned short) + ARRAY_READERS(int) + ARRAY_READERS(unsigned int) + ARRAY_READERS(unsigned long int) + ARRAY_READERS(double) ARRAY_READERS(float) #undef ARRAY_READERS +#define NDFIELD_READERS(D, N) \ + template void ReadNDField(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + ndfield_t&, \ + const adios2::Box&); + NDFIELD_READERS(Dim::_1D, 3) + NDFIELD_READERS(Dim::_1D, 6) + NDFIELD_READERS(Dim::_2D, 3) + NDFIELD_READERS(Dim::_2D, 6) + NDFIELD_READERS(Dim::_3D, 3) + NDFIELD_READERS(Dim::_3D, 6) +#undef NDFIELD_READERS + } // namespace out diff --git a/src/output/utils/readers.h b/src/output/utils/readers.h index b4fdc8ab0..c62489282 100644 --- a/src/output/utils/readers.h +++ b/src/output/utils/readers.h @@ -6,6 +6,7 @@ * - out::ReadVariable<> -> void * - out::Read1DArray<> -> void * - out::Read2DArray<> -> void + * - out::ReadNDField<> -> void * @cpp: * - readers.cpp * @namespaces: @@ -43,6 +44,13 @@ namespace out { std::size_t, std::size_t); + template + void ReadNDField(adios2::IO&, + adios2::Engine&, + const std::string&, + ndfield_t&, + const adios2::Box&); + } // namespace out #endif // OUTPUT_UTILS_READERS_H diff --git a/src/output/utils/writers.cpp b/src/output/utils/writers.cpp index c7d1bbd74..b16ca9e8d 100644 --- a/src/output/utils/writers.cpp +++ b/src/output/utils/writers.cpp @@ -15,7 +15,7 @@ namespace out { const T& data, std::size_t global_size, std::size_t local_offset) { - auto var = io.InquireVariable(name); + auto var = io.InquireVariable(name); var.SetShape({ global_size }); var.SetSelection(adios2::Box({ local_offset }, { 1 })); writer.Put(var, &data); @@ -27,8 +27,8 @@ namespace out { const std::string& name, const array_t& data, std::size_t local_size, - std::size_t local_offset, - std::size_t global_size) { + std::size_t global_size, + std::size_t local_offset) { const auto slice = range_tuple_t(0, local_size); auto var = io.InquireVariable(name); var.SetShape({ global_size }); @@ -47,8 +47,8 @@ namespace out { const array_t& data, unsigned short dim2_size, std::size_t local_size, - std::size_t local_offset, - std::size_t global_size) { + std::size_t global_size, + std::size_t local_offset) { const auto slice = range_tuple_t(0, local_size); auto var = io.InquireVariable(name); @@ -62,6 +62,16 @@ namespace out { writer.Put(var, data_sub.data(), adios2::Mode::Sync); } + template + void WriteNDField(adios2::IO& io, + adios2::Engine& writer, + const std::string& name, + const ndfield_t& data) { + auto data_h = Kokkos::create_mirror_view(data); + Kokkos::deep_copy(data_h, data); + writer.Put(io.InquireVariable(name), data_h.data(), adios2::Mode::Sync); + } + #define ARRAY_WRITERS(T) \ template void WriteVariable(adios2::IO&, \ adios2::Engine&, \ @@ -83,12 +93,28 @@ namespace out { unsigned short, \ std::size_t, \ std::size_t, \ - std::size_t); \ - ARRAY_WRITERS(int) \ - ARRAY_WRITERS(unsigned int) \ - ARRAY_WRITERS(unsigned long int) \ - ARRAY_WRITERS(double) \ + std::size_t); + + ARRAY_WRITERS(short) + ARRAY_WRITERS(unsigned short) + ARRAY_WRITERS(int) + ARRAY_WRITERS(unsigned int) + ARRAY_WRITERS(unsigned long int) + ARRAY_WRITERS(double) ARRAY_WRITERS(float) #undef ARRAY_WRITERS +#define NDFIELD_WRITERS(D, N) \ + template void WriteNDField(adios2::IO&, \ + adios2::Engine&, \ + const std::string&, \ + const ndfield_t&); + NDFIELD_WRITERS(Dim::_1D, 3) + NDFIELD_WRITERS(Dim::_1D, 6) + NDFIELD_WRITERS(Dim::_2D, 3) + NDFIELD_WRITERS(Dim::_2D, 6) + NDFIELD_WRITERS(Dim::_3D, 3) + NDFIELD_WRITERS(Dim::_3D, 6) +#undef NDFIELD_WRITERS + } // namespace out diff --git a/src/output/utils/writers.h b/src/output/utils/writers.h index 3e4bd4cb8..58fac2bf6 100644 --- a/src/output/utils/writers.h +++ b/src/output/utils/writers.h @@ -6,6 +6,7 @@ * - out::WriteVariable<> -> void * - out::Write1DArray<> -> void * - out::Write2DArray<> -> void + * - out::WriteNDField<> -> void * @cpp: * - writers.cpp * @namespaces: @@ -50,6 +51,12 @@ namespace out { std::size_t, std::size_t); + template + void WriteNDField(adios2::IO&, + adios2::Engine&, + const std::string&, + const ndfield_t&); + } // namespace out #endif // OUTPUT_UTILS_WRITERS_H From ce4e96f9117fe818188d549550f22f7c39b3a93b Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 21 Oct 2025 16:09:07 -0400 Subject: [PATCH 06/15] minor (rm ntt_checkpoint) --- cmake/benchmark.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/benchmark.cmake b/cmake/benchmark.cmake index 39b075716..fdd8438ea 100644 --- a/cmake/benchmark.cmake +++ b/cmake/benchmark.cmake @@ -20,7 +20,7 @@ add_executable(${exec} ${src}) set(libs ntt_global ntt_metrics ntt_kernels ntt_archetypes ntt_framework) if(${output}) - list(APPEND libs ntt_output ntt_checkpoint) + list(APPEND libs ntt_output) endif() add_dependencies(${exec} ${libs}) target_link_libraries(${exec} PRIVATE ${libs} stdc++fs) From 662b7f7677d684386c7e1af926977f03ccd3c942 Mon Sep 17 00:00:00 2001 From: hayk Date: Fri, 24 Oct 2025 15:49:36 -0400 Subject: [PATCH 07/15] experimental removed (WIP) --- src/archetypes/energy_dist.h | 309 ++++++----------- src/archetypes/particle_injector.h | 514 ++++++++++++++--------------- src/archetypes/utils.h | 39 +-- src/kernels/injectors.hpp | 387 ++++++++++------------ 4 files changed, 541 insertions(+), 708 deletions(-) diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h index bcb88fbc3..f8ef1175a 100644 --- a/src/archetypes/energy_dist.h +++ b/src/archetypes/energy_dist.h @@ -208,228 +208,125 @@ namespace arch { struct Maxwellian : public EnergyDistribution { using EnergyDistribution::metric; - Maxwellian(const M& metric, - random_number_pool_t& pool, - real_t temperature, - real_t boost_vel = ZERO, - in boost_direction = in::x1, - bool zero_current = true) + Maxwellian(const M& metric, + random_number_pool_t& pool, + real_t temperature, + const std::vector& drift_four_vel = { ZERO, ZERO, ZERO }) : EnergyDistribution { metric } , pool { pool } - , temperature { temperature } - , boost_velocity { boost_vel } - , boost_direction { boost_direction } - , zero_current { zero_current } { + , temperature { temperature } { + raise::ErrorIf(drift_four_vel.size() != 3, + "Maxwellian: Drift velocity must be a 3D vector", + HERE); raise::ErrorIf(temperature < ZERO, "Maxwellian: Temperature must be non-negative", HERE); - raise::ErrorIf( - (not cmp::AlmostZero_host(boost_vel, ZERO)) && (M::CoordType != Coord::Cart), - "Maxwellian: Boosting is only supported in Cartesian coordinates", - HERE); - } - - Inline void operator()(const coord_t&, - vec_t& v, - spidx_t sp = 0) const { - SampleFromMaxwellian(v, - pool, - temperature, - boost_velocity, - boost_direction, - not zero_current and - sp % 2 == 0); - } - - private: - random_number_pool_t pool; - - const real_t temperature; - const real_t boost_velocity; - const in boost_direction; - const bool zero_current; - }; - - template - struct TwoTemperatureMaxwellian : public EnergyDistribution { - using EnergyDistribution::metric; - - TwoTemperatureMaxwellian(const M& metric, - random_number_pool_t& pool, - const std::pair& temperatures, - const std::pair& species, - real_t boost_vel = ZERO, - in boost_direction = in::x1, - bool zero_current = true) - : EnergyDistribution { metric } - , pool { pool } - , temperature_1 { temperatures.first } - , temperature_2 { temperatures.second } - , sp_1 { species.first } - , sp_2 { species.second } - , boost_velocity { boost_vel } - , boost_direction { boost_direction } - , zero_current { zero_current } { - raise::ErrorIf( - (temperature_1 < ZERO) or (temperature_2 < ZERO), - "TwoTemperatureMaxwellian: Temperature must be non-negative", - HERE); - raise::ErrorIf((not cmp::AlmostZero(boost_vel, ZERO)) && - (M::CoordType != Coord::Cart), - "TwoTemperatureMaxwellian: Boosting is only supported in " - "Cartesian coordinates", - HERE); - } - - Inline void operator()(const coord_t&, - vec_t& v, - spidx_t sp = 0) const { - SampleFromMaxwellian( - v, - pool, - (sp == sp_1) ? temperature_1 : temperature_2, - boost_velocity, - boost_direction, - not zero_current and sp == sp_1); - } - - private: - random_number_pool_t pool; - - const real_t temperature_1, temperature_2; - const spidx_t sp_1, sp_2; - const real_t boost_velocity; - const in boost_direction; - const bool zero_current; - }; - - namespace experimental { - - template - struct Maxwellian : public EnergyDistribution { - using EnergyDistribution::metric; - - Maxwellian(const M& metric, - random_number_pool_t& pool, - real_t temperature, - const std::vector& drift_four_vel = { ZERO, ZERO, ZERO }) - : EnergyDistribution { metric } - , pool { pool } - , temperature { temperature } { - raise::ErrorIf(drift_four_vel.size() != 3, - "Maxwellian: Drift velocity must be a 3D vector", - HERE); - raise::ErrorIf(temperature < ZERO, - "Maxwellian: Temperature must be non-negative", - HERE); - if constexpr (M::CoordType == Coord::Cart) { - drift_4vel = NORM(drift_four_vel[0], drift_four_vel[1], drift_four_vel[2]); - if (cmp::AlmostZero_host(drift_4vel)) { - drift_dir = 0; - } else { - drift_3vel = drift_4vel / math::sqrt(ONE + SQR(drift_4vel)); - drift_dir_x1 = drift_four_vel[0] / drift_4vel; - drift_dir_x2 = drift_four_vel[1] / drift_4vel; - drift_dir_x3 = drift_four_vel[2] / drift_4vel; - - // assume drift is in an arbitrary direction - drift_dir = 4; - // check whether drift is in one of principal directions - for (auto d { 0u }; d < 3u; ++d) { - const auto dprev = (d + 2) % 3; - const auto dnext = (d + 1) % 3; - if (cmp::AlmostZero_host(drift_four_vel[dprev]) and - cmp::AlmostZero_host(drift_four_vel[dnext])) { - drift_dir = SIGN(drift_four_vel[d]) * (d + 1); - break; - } + if constexpr (M::CoordType == Coord::Cart) { + drift_4vel = NORM(drift_four_vel[0], drift_four_vel[1], drift_four_vel[2]); + if (cmp::AlmostZero_host(drift_4vel)) { + drift_dir = 0; + } else { + drift_3vel = drift_4vel / math::sqrt(ONE + SQR(drift_4vel)); + drift_dir_x1 = drift_four_vel[0] / drift_4vel; + drift_dir_x2 = drift_four_vel[1] / drift_4vel; + drift_dir_x3 = drift_four_vel[2] / drift_4vel; + + // assume drift is in an arbitrary direction + drift_dir = 4; + // check whether drift is in one of principal directions + for (auto d { 0u }; d < 3u; ++d) { + const auto dprev = (d + 2) % 3; + const auto dnext = (d + 1) % 3; + if (cmp::AlmostZero_host(drift_four_vel[dprev]) and + cmp::AlmostZero_host(drift_four_vel[dnext])) { + drift_dir = SIGN(drift_four_vel[d]) * (d + 1); + break; } } - raise::ErrorIf(drift_dir > 3 and drift_dir != 4, - "Maxwellian: Incorrect drift direction", - HERE); - raise::ErrorIf( - drift_dir != 0 and (M::CoordType != Coord::Cart), - "Maxwellian: Boosting is only supported in Cartesian coordinates", - HERE); } + raise::ErrorIf(drift_dir > 3 and drift_dir != 4, + "Maxwellian: Incorrect drift direction", + HERE); + raise::ErrorIf( + drift_dir != 0 and (M::CoordType != Coord::Cart), + "Maxwellian: Boosting is only supported in Cartesian coordinates", + HERE); } + } - Inline void operator()(const coord_t& x_Code, - vec_t& v, - spidx_t = 0) const { - if (cmp::AlmostZero(temperature)) { - v[0] = ZERO; - v[1] = ZERO; - v[2] = ZERO; - } else { - JuttnerSinge(v, temperature, pool); - } - // @note: boost only when using cartesian coordinates - if constexpr (M::CoordType == Coord::Cart) { - if (drift_dir != 0) { - // Boost an isotropic Maxwellian with a drift velocity using - // flipping method https://arxiv.org/pdf/1504.03910.pdf - // 1. apply drift in X1 direction - const auto gamma { U2GAMMA(v[0], v[1], v[2]) }; - auto rand_gen = pool.get_state(); - if (-drift_3vel * v[0] > gamma * Random(rand_gen)) { - v[0] = -v[0]; - } - pool.free_state(rand_gen); - v[0] = math::sqrt(ONE + SQR(drift_4vel)) * (v[0] + drift_3vel * gamma); - // 2. rotate to desired orientation - if (drift_dir == -1) { - v[0] = -v[0]; - } else if (drift_dir == 2 || drift_dir == -2) { - const auto tmp = v[1]; - v[1] = drift_dir > 0 ? v[0] : -v[0]; - v[0] = tmp; - } else if (drift_dir == 3 || drift_dir == -3) { - const auto tmp = v[2]; - v[2] = drift_dir > 0 ? v[0] : -v[0]; - v[0] = tmp; - } else if (drift_dir == 4) { - vec_t v_old; - v_old[0] = v[0]; - v_old[1] = v[1]; - v_old[2] = v[2]; - - v[0] = v_old[0] * drift_dir_x1 - v_old[1] * drift_dir_x2 - - v_old[2] * drift_dir_x3; - v[1] = (v_old[0] * drift_dir_x2 * (drift_dir_x1 + ONE) + - v_old[1] * - (SQR(drift_dir_x1) + drift_dir_x1 + SQR(drift_dir_x3)) - - v_old[2] * drift_dir_x2 * drift_dir_x3) / - (drift_dir_x1 + ONE); - v[2] = (v_old[0] * drift_dir_x3 * (drift_dir_x1 + ONE) - - v_old[1] * drift_dir_x2 * drift_dir_x3 - - v_old[2] * (-drift_dir_x1 + SQR(drift_dir_x3) - ONE)) / - (drift_dir_x1 + ONE); - } + Inline void operator()(const coord_t& x_Code, + vec_t& v, + spidx_t = 0) const { + if (cmp::AlmostZero(temperature)) { + v[0] = ZERO; + v[1] = ZERO; + v[2] = ZERO; + } else { + JuttnerSinge(v, temperature, pool); + } + // @note: boost only when using cartesian coordinates + if constexpr (M::CoordType == Coord::Cart) { + if (drift_dir != 0) { + // Boost an isotropic Maxwellian with a drift velocity using + // flipping method https://arxiv.org/pdf/1504.03910.pdf + // 1. apply drift in X1 direction + const auto gamma { U2GAMMA(v[0], v[1], v[2]) }; + auto rand_gen = pool.get_state(); + if (-drift_3vel * v[0] > gamma * Random(rand_gen)) { + v[0] = -v[0]; + } + pool.free_state(rand_gen); + v[0] = math::sqrt(ONE + SQR(drift_4vel)) * (v[0] + drift_3vel * gamma); + // 2. rotate to desired orientation + if (drift_dir == -1) { + v[0] = -v[0]; + } else if (drift_dir == 2 || drift_dir == -2) { + const auto tmp = v[1]; + v[1] = drift_dir > 0 ? v[0] : -v[0]; + v[0] = tmp; + } else if (drift_dir == 3 || drift_dir == -3) { + const auto tmp = v[2]; + v[2] = drift_dir > 0 ? v[0] : -v[0]; + v[0] = tmp; + } else if (drift_dir == 4) { + vec_t v_old; + v_old[0] = v[0]; + v_old[1] = v[1]; + v_old[2] = v[2]; + + v[0] = v_old[0] * drift_dir_x1 - v_old[1] * drift_dir_x2 - + v_old[2] * drift_dir_x3; + v[1] = (v_old[0] * drift_dir_x2 * (drift_dir_x1 + ONE) + + v_old[1] * + (SQR(drift_dir_x1) + drift_dir_x1 + SQR(drift_dir_x3)) - + v_old[2] * drift_dir_x2 * drift_dir_x3) / + (drift_dir_x1 + ONE); + v[2] = (v_old[0] * drift_dir_x3 * (drift_dir_x1 + ONE) - + v_old[1] * drift_dir_x2 * drift_dir_x3 - + v_old[2] * (-drift_dir_x1 + SQR(drift_dir_x3) - ONE)) / + (drift_dir_x1 + ONE); } } } + } - private: - random_number_pool_t pool; - - const real_t temperature; - - real_t drift_3vel { ZERO }, drift_4vel { ZERO }; - // components of the unit vector in the direction of the drift - real_t drift_dir_x1 { ZERO }, drift_dir_x2 { ZERO }, drift_dir_x3 { ZERO }; + private: + random_number_pool_t pool; - // values of boost_dir: - // 4 -> arbitrary direction - // 0 -> no drift - // +/- 1 -> +/- x1 - // +/- 2 -> +/- x2 - // +/- 3 -> +/- x3 - short drift_dir { 0 }; - }; + const real_t temperature; - } // namespace experimental + real_t drift_3vel { ZERO }, drift_4vel { ZERO }; + // components of the unit vector in the direction of the drift + real_t drift_dir_x1 { ZERO }, drift_dir_x2 { ZERO }, drift_dir_x3 { ZERO }; + + // values of boost_dir: + // 4 -> arbitrary direction + // 0 -> no drift + // +/- 1 -> +/- x1 + // +/- 2 -> +/- x2 + // +/- 3 -> +/- x3 + short drift_dir { 0 }; + }; } // namespace arch diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 6313031d1..f92f6defa 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -28,7 +28,6 @@ #include "framework/domain/metadomain.h" #include "kernels/injectors.hpp" -#include "kernels/particle_moments.hpp" #include "kernels/utils.hpp" #include @@ -117,149 +116,146 @@ namespace arch { } }; - template class ED> - struct UniformInjector : BaseInjector { - using energy_dist_t = ED; - static_assert(M::is_metric, "M must be a metric class"); - static_assert(energy_dist_t::is_energy_dist, - "E must be an energy distribution class"); - static constexpr bool is_uniform_injector { true }; - static constexpr Dimension D { M::Dim }; - static constexpr Coord C { M::CoordType }; - - const energy_dist_t energy_dist; - const std::pair species; - - UniformInjector(const energy_dist_t& energy_dist, - const std::pair& species) - : energy_dist { energy_dist } - , species { species } {} - - ~UniformInjector() = default; - }; - - template class ED> - struct KeepConstantInjector : UniformInjector { - using energy_dist_t = ED; - using UniformInjector::D; - using UniformInjector::C; - - const idx_t density_buff_idx; - boundaries_t probe_box; - - KeepConstantInjector(const energy_dist_t& energy_dist, - const std::pair& species, - idx_t density_buff_idx, - boundaries_t box = {}) - : UniformInjector { energy_dist, species } - , density_buff_idx { density_buff_idx } { - for (auto d { 0u }; d < M::Dim; ++d) { - if (d < box.size()) { - probe_box.push_back({ box[d].first, box[d].second }); - } else { - probe_box.push_back(Range::All); - } - } - } - - ~KeepConstantInjector() = default; - - auto ComputeAvgDensity(const SimulationParams& params, - const Domain& domain) const -> real_t { - const auto result = this->DeduceRegion(domain, probe_box); - const auto should_probe = std::get<0>(result); - if (not should_probe) { - return ZERO; - } - const auto xi_min_arr = std::get<1>(result); - const auto xi_max_arr = std::get<2>(result); - - tuple_t i_min { 0 }; - tuple_t i_max { 0 }; - - auto xi_min_h = Kokkos::create_mirror_view(xi_min_arr); - auto xi_max_h = Kokkos::create_mirror_view(xi_max_arr); - Kokkos::deep_copy(xi_min_h, xi_min_arr); - Kokkos::deep_copy(xi_max_h, xi_max_arr); - - ncells_t num_cells = 1u; - for (auto d { 0u }; d < M::Dim; ++d) { - i_min[d] = std::floor(xi_min_h(d)) + N_GHOSTS; - i_max[d] = std::ceil(xi_max_h(d)) + N_GHOSTS; - num_cells *= (i_max[d] - i_min[d]); - } - - real_t dens { ZERO }; - if (should_probe) { - Kokkos::parallel_reduce( - "AvgDensity", - CreateRangePolicy(i_min, i_max), - kernel::ComputeSum_kernel(domain.fields.buff, density_buff_idx), - dens); - } -#if defined(MPI_ENABLED) - real_t tot_dens { ZERO }; - ncells_t tot_num_cells { 0 }; - MPI_Allreduce(&dens, &tot_dens, 1, mpi::get_type(), MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&num_cells, - &tot_num_cells, - 1, - mpi::get_type(), - MPI_SUM, - MPI_COMM_WORLD); - dens = tot_dens; - num_cells = tot_num_cells; -#endif - if (num_cells > 0) { - return dens / (real_t)(num_cells); - } else { - return ZERO; - } - } - - auto ComputeNumInject(const SimulationParams& params, - const Domain& domain, - real_t number_density, - const boundaries_t& box) const - -> std::tuple, array_t> override { - const auto computed_avg_density = ComputeAvgDensity(params, domain); - - const auto result = this->DeduceRegion(domain, box); - if (not std::get<0>(result)) { - return { false, (npart_t)0, array_t {}, array_t {} }; - } - - const auto xi_min = std::get<1>(result); - const auto xi_max = std::get<2>(result); - auto xi_min_h = Kokkos::create_mirror_view(xi_min); - auto xi_max_h = Kokkos::create_mirror_view(xi_max); - Kokkos::deep_copy(xi_min_h, xi_min); - Kokkos::deep_copy(xi_max_h, xi_max); - - long double num_cells { 1.0 }; - for (auto d { 0u }; d < M::Dim; ++d) { - num_cells *= static_cast(xi_max_h(d)) - - static_cast(xi_min_h(d)); - } - - const auto ppc0 = params.template get("particles.ppc0"); - npart_t nparticles { 0u }; - if (number_density > computed_avg_density) { - nparticles = static_cast( - (long double)(ppc0 * (number_density - computed_avg_density) * 0.5) * - num_cells); - } - - return { nparticles != 0u, nparticles, xi_min, xi_max }; - } - }; + // template + // class ED> struct UniformInjector : BaseInjector { + // using energy_dist_t = ED; + // static_assert(M::is_metric, "M must be a metric class"); + // static_assert(energy_dist_t::is_energy_dist, + // "E must be an energy distribution class"); + // static constexpr bool is_uniform_injector { true }; + // static constexpr Dimension D { M::Dim }; + // static constexpr Coord C { M::CoordType }; + // + // const energy_dist_t energy_dist; + // const std::pair species; + // + // UniformInjector(const energy_dist_t& energy_dist, + // const std::pair& species) + // : energy_dist { energy_dist } + // , species { species } {} + // + // ~UniformInjector() = default; + // }; + + // template class ED> + // struct KeepConstantInjector : UniformInjector { + // using energy_dist_t = ED; + // using UniformInjector::D; + // using UniformInjector::C; + // + // const idx_t density_buff_idx; + // boundaries_t probe_box; + // + // KeepConstantInjector(const energy_dist_t& energy_dist, + // const std::pair& species, + // idx_t density_buff_idx, boundaries_t box = {}) + // : UniformInjector { energy_dist, species } + // , density_buff_idx { density_buff_idx } { + // for (auto d { 0u }; d < M::Dim; ++d) { + // if (d < box.size()) { + // probe_box.push_back({ box[d].first, box[d].second }); + // } else { + // probe_box.push_back(Range::All); + // } + // } + // } + // + // ~KeepConstantInjector() = default; + // + // auto ComputeAvgDensity(const SimulationParams& params, + // const Domain& domain) const -> real_t { + // const auto result = this->DeduceRegion(domain, probe_box); + // const auto should_probe = std::get<0>(result); + // if (not should_probe) { + // return ZERO; + // } + // const auto xi_min_arr = std::get<1>(result); + // const auto xi_max_arr = std::get<2>(result); + // + // tuple_t i_min { 0 }; + // tuple_t i_max { 0 }; + // + // auto xi_min_h = Kokkos::create_mirror_view(xi_min_arr); + // auto xi_max_h = Kokkos::create_mirror_view(xi_max_arr); + // Kokkos::deep_copy(xi_min_h, xi_min_arr); + // Kokkos::deep_copy(xi_max_h, xi_max_arr); + // + // ncells_t num_cells = 1u; + // for (auto d { 0u }; d < M::Dim; ++d) { + // i_min[d] = std::floor(xi_min_h(d)) + N_GHOSTS; + // i_max[d] = std::ceil(xi_max_h(d)) + N_GHOSTS; + // num_cells *= (i_max[d] - i_min[d]); + // } + // + // real_t dens { ZERO }; + // if (should_probe) { + // Kokkos::parallel_reduce( + // "AvgDensity", + // CreateRangePolicy(i_min, i_max), + // kernel::ComputeSum_kernel(domain.fields.buff, density_buff_idx), + // dens); + // } + // #if defined(MPI_ENABLED) + // real_t tot_dens { ZERO }; + // ncells_t tot_num_cells { 0 }; + // MPI_Allreduce(&dens, &tot_dens, 1, mpi::get_type(), MPI_SUM, MPI_COMM_WORLD); + // MPI_Allreduce(&num_cells, + // &tot_num_cells, + // 1, + // mpi::get_type(), + // MPI_SUM, + // MPI_COMM_WORLD); + // dens = tot_dens; + // num_cells = tot_num_cells; + // #endif + // if (num_cells > 0) { + // return dens / (real_t)(num_cells); + // } else { + // return ZERO; + // } + // } + // + // auto ComputeNumInject(const SimulationParams& params, + // const Domain& domain, + // real_t number_density, + // const boundaries_t& box) const + // -> std::tuple, array_t> override { + // const auto computed_avg_density = ComputeAvgDensity(params, domain); + // + // const auto result = this->DeduceRegion(domain, box); + // if (not std::get<0>(result)) { + // return { false, (npart_t)0, array_t {}, array_t {} }; + // } + // + // const auto xi_min = std::get<1>(result); + // const auto xi_max = std::get<2>(result); + // auto xi_min_h = Kokkos::create_mirror_view(xi_min); + // auto xi_max_h = Kokkos::create_mirror_view(xi_max); + // Kokkos::deep_copy(xi_min_h, xi_min); + // Kokkos::deep_copy(xi_max_h, xi_max); + // + // long double num_cells { 1.0 }; + // for (auto d { 0u }; d < M::Dim; ++d) { + // num_cells *= static_cast(xi_max_h(d)) - + // static_cast(xi_min_h(d)); + // } + // + // const auto ppc0 = params.template get("particles.ppc0"); + // npart_t nparticles { 0u }; + // if (number_density > computed_avg_density) { + // nparticles = static_cast( + // (long double)(ppc0 * (number_density - computed_avg_density) * 0.5) * + // num_cells); + // } + // + // return { nparticles != 0u, nparticles, xi_min, xi_max }; + // } + // }; template - class ED, - template - class SD> + template class ED, + template class SD> struct NonUniformInjector { using energy_dist_t = ED; using spatial_dist_t = SD; @@ -426,6 +422,117 @@ namespace arch { ~MovingInjector() = default; }; + // /** + // * @brief Injects uniform number density of particles everywhere in the domain + // * @param domain Domain object + // * @param injector Uniform injector object + // * @param number_density Total number density (in units of n0) + // * @param use_weights Use weights + // * @param box Region to inject the particles in global coords + // * @tparam S Simulation engine type + // * @tparam M Metric type + // * @tparam I Injector type + // */ + // template + // inline void InjectUniform(const SimulationParams& params, + // Domain& domain, + // const I& injector, + // real_t number_density, + // bool use_weights = false, + // const boundaries_t& box = {}) { + // static_assert(M::is_metric, "M must be a metric class"); + // static_assert(I::is_uniform_injector, "I must be a uniform injector class"); + // raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights), + // "Weights must be used for non-Cartesian coordinates", + // HERE); + // raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights, + // "Weights should not be used for Cartesian coordinates", + // HERE); + // raise::ErrorIf(params.template get("particles.use_weights") != use_weights, + // "Weights must be enabled from the input file to use them in " + // "the injector", + // HERE); + // if (domain.species[injector.species.first - 1].charge() + + // domain.species[injector.species.second - 1].charge() != + // 0.0f) { + // raise::Warning("Total charge of the injected species is non-zero", HERE); + // } + // + // { + // boundaries_t nonempty_box; + // for (auto d { 0u }; d < M::Dim; ++d) { + // if (d < box.size()) { + // nonempty_box.push_back({ box[d].first, box[d].second }); + // } else { + // nonempty_box.push_back(Range::All); + // } + // } + // const auto result = injector.ComputeNumInject(params, + // domain, + // number_density, + // nonempty_box); + // if (not std::get<0>(result)) { + // return; + // } + // const auto nparticles = std::get<1>(result); + // const auto xi_min = std::get<2>(result); + // const auto xi_max = std::get<3>(result); + // + // Kokkos::parallel_for( + // "InjectUniform", + // nparticles, + // kernel::UniformInjector_kernel( + // injector.species.first, + // injector.species.second, + // domain.species[injector.species.first - 1], + // domain.species[injector.species.second - 1], + // domain.species[injector.species.first - 1].npart(), + // domain.species[injector.species.second - 1].npart(), + // domain.mesh.metric, + // xi_min, + // xi_max, + // injector.energy_dist, + // ONE / params.template get("scales.V0"), + // domain.random_pool)); + // domain.species[injector.species.first - 1].set_npart( + // domain.species[injector.species.first - 1].npart() + nparticles); + // domain.species[injector.species.second - 1].set_npart( + // domain.species[injector.species.second - 1].npart() + nparticles); + // } + // } + // + // namespace experimental { + + template class ED1, + template class ED2> + struct UniformInjector : BaseInjector { + using energy_dist_1_t = ED1; + using energy_dist_2_t = ED2; + static_assert(M::is_metric, "M must be a metric class"); + static_assert(energy_dist_1_t::is_energy_dist, + "ED1 must be an energy distribution class"); + static_assert(energy_dist_2_t::is_energy_dist, + "ED2 must be an energy distribution class"); + static constexpr bool is_uniform_injector { true }; + static constexpr Dimension D { M::Dim }; + static constexpr Coord C { M::CoordType }; + + const energy_dist_1_t energy_dist_1; + const energy_dist_2_t energy_dist_2; + const std::pair species; + + UniformInjector(const energy_dist_1_t& energy_dist_1, + const energy_dist_2_t& energy_dist_2, + const std::pair& species) + : energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } + , species { species } {} + + ~UniformInjector() = default; + }; + /** * @brief Injects uniform number density of particles everywhere in the domain * @param domain Domain object @@ -485,17 +592,20 @@ namespace arch { Kokkos::parallel_for( "InjectUniform", nparticles, - kernel::UniformInjector_kernel( + kernel::UniformInjector_kernel( injector.species.first, injector.species.second, domain.species[injector.species.first - 1], domain.species[injector.species.second - 1], + nparticles, + domain.index(), domain.species[injector.species.first - 1].npart(), domain.species[injector.species.second - 1].npart(), domain.mesh.metric, xi_min, xi_max, - injector.energy_dist, + injector.energy_dist_1, + injector.energy_dist_2, ONE / params.template get("scales.V0"), domain.random_pool)); domain.species[injector.species.first - 1].set_npart( @@ -505,123 +615,7 @@ namespace arch { } } - namespace experimental { - - template - class ED1, - template - class ED2> - struct UniformInjector : BaseInjector { - using energy_dist_1_t = ED1; - using energy_dist_2_t = ED2; - static_assert(M::is_metric, "M must be a metric class"); - static_assert(energy_dist_1_t::is_energy_dist, - "ED1 must be an energy distribution class"); - static_assert(energy_dist_2_t::is_energy_dist, - "ED2 must be an energy distribution class"); - static constexpr bool is_uniform_injector { true }; - static constexpr Dimension D { M::Dim }; - static constexpr Coord C { M::CoordType }; - - const energy_dist_1_t energy_dist_1; - const energy_dist_2_t energy_dist_2; - const std::pair species; - - UniformInjector(const energy_dist_1_t& energy_dist_1, - const energy_dist_2_t& energy_dist_2, - const std::pair& species) - : energy_dist_1 { energy_dist_1 } - , energy_dist_2 { energy_dist_2 } - , species { species } {} - - ~UniformInjector() = default; - }; - - /** - * @brief Injects uniform number density of particles everywhere in the domain - * @param domain Domain object - * @param injector Uniform injector object - * @param number_density Total number density (in units of n0) - * @param use_weights Use weights - * @param box Region to inject the particles in global coords - * @tparam S Simulation engine type - * @tparam M Metric type - * @tparam I Injector type - */ - template - inline void InjectUniform(const SimulationParams& params, - Domain& domain, - const I& injector, - real_t number_density, - bool use_weights = false, - const boundaries_t& box = {}) { - static_assert(M::is_metric, "M must be a metric class"); - static_assert(I::is_uniform_injector, "I must be a uniform injector class"); - raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights), - "Weights must be used for non-Cartesian coordinates", - HERE); - raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights, - "Weights should not be used for Cartesian coordinates", - HERE); - raise::ErrorIf( - params.template get("particles.use_weights") != use_weights, - "Weights must be enabled from the input file to use them in " - "the injector", - HERE); - if (domain.species[injector.species.first - 1].charge() + - domain.species[injector.species.second - 1].charge() != - 0.0f) { - raise::Warning("Total charge of the injected species is non-zero", HERE); - } - - { - boundaries_t nonempty_box; - for (auto d { 0u }; d < M::Dim; ++d) { - if (d < box.size()) { - nonempty_box.push_back({ box[d].first, box[d].second }); - } else { - nonempty_box.push_back(Range::All); - } - } - const auto result = injector.ComputeNumInject(params, - domain, - number_density, - nonempty_box); - if (not std::get<0>(result)) { - return; - } - const auto nparticles = std::get<1>(result); - const auto xi_min = std::get<2>(result); - const auto xi_max = std::get<3>(result); - - Kokkos::parallel_for( - "InjectUniform", - nparticles, - kernel::experimental:: - UniformInjector_kernel( - injector.species.first, - injector.species.second, - domain.species[injector.species.first - 1], - domain.species[injector.species.second - 1], - domain.species[injector.species.first - 1].npart(), - domain.species[injector.species.second - 1].npart(), - domain.mesh.metric, - xi_min, - xi_max, - injector.energy_dist_1, - injector.energy_dist_2, - ONE / params.template get("scales.V0"), - domain.random_pool)); - domain.species[injector.species.first - 1].set_npart( - domain.species[injector.species.first - 1].npart() + nparticles); - domain.species[injector.species.second - 1].set_npart( - domain.species[injector.species.second - 1].npart() + nparticles); - } - } - - } // namespace experimental + // } // namespace experimental /** * @brief Injects particles from a globally-defined map diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h index d3d6ae475..7a5296771 100644 --- a/src/archetypes/utils.h +++ b/src/archetypes/utils.h @@ -51,29 +51,26 @@ namespace arch { const auto temperature_1 = temperatures.first / mass_1; const auto temperature_2 = temperatures.second / mass_2; - const auto maxwellian_1 = arch::experimental::Maxwellian( - domain.mesh.metric, - domain.random_pool, - temperature_1, - drift_four_vels.first); - const auto maxwellian_2 = arch::experimental::Maxwellian( - domain.mesh.metric, - domain.random_pool, - temperature_2, - drift_four_vels.second); + const auto maxwellian_1 = arch::Maxwellian(domain.mesh.metric, + domain.random_pool, + temperature_1, + drift_four_vels.first); + const auto maxwellian_2 = arch::Maxwellian(domain.mesh.metric, + domain.random_pool, + temperature_2, + drift_four_vels.second); - const auto injector = arch::experimental:: - UniformInjector( - maxwellian_1, - maxwellian_2, - species); + const auto injector = arch::UniformInjector( + maxwellian_1, + maxwellian_2, + species); - arch::experimental::InjectUniform(params, - domain, - injector, - tot_number_density, - use_weights, - box); + arch::InjectUniform(params, + domain, + injector, + tot_number_density, + use_weights, + box); } /** diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index c321d18f8..69642e2fd 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -25,9 +25,60 @@ namespace kernel { using namespace ntt; - template + template + Inline void InjectParticle(const npart_t& p, + const array_t& i1_arr, + const array_t& i2_arr, + const array_t& i3_arr, + const array_t& dx1_arr, + const array_t& dx2_arr, + const array_t& dx3_arr, + const array_t& ux1_arr, + const array_t& ux2_arr, + const array_t& ux3_arr, + const array_t& phi_arr, + const array_t& weight_arr, + const array_t& tag_arr, + const array_t& pld_i_arr, + const tuple_t& xi_Cd, + const tuple_t& dxi_Cd, + const vec_t& v_Cd, + const real_t& weight = ONE, + const real_t& phi = ZERO, + const npart_t& domain_idx = 0u, + const npart_t& species_cntr = 0u) { + if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { + i1_arr(p) = xi_Cd[0]; + dx1_arr(p) = dxi_Cd[0]; + } + if constexpr (D == Dim::_2D or D == Dim::_3D) { + i2_arr(p) = xi_Cd[1]; + dx2_arr(p) = dxi_Cd[1]; + } + if constexpr (D == Dim::_3D) { + i3_arr(p) = xi_Cd[2]; + dx3_arr(p) = dxi_Cd[2]; + } + if constexpr (D == Dim::_2D and C != Coord::Cart) { + phi_arr(p) = phi; + } + ux1_arr(p) = v_Cd[0]; + ux2_arr(p) = v_Cd[1]; + ux3_arr(p) = v_Cd[2]; + tag_arr(p) = ParticleTag::alive; + weight_arr(p) = weight; + if constexpr (T) { + pld_i_arr(p, pldi::spcCtr) = species_cntr; +#if defined(MPI_ENABLED) + pld_i_arr(p, pldi::domIdx) = domain_idx; +#endif + } + } + + template struct UniformInjector_kernel { - static_assert(ED::is_energy_dist, "ED must be an energy distribution class"); + static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class"); static_assert(M::is_metric, "M must be a metric class"); const spidx_t spidx1, spidx2; @@ -38,6 +89,7 @@ namespace kernel { array_t phis_1; array_t weights_1; array_t tags_1; + array_t pldis_1; array_t i1s_2, i2s_2, i3s_2; array_t dx1s_2, dx2s_2, dx3s_2; @@ -45,11 +97,15 @@ namespace kernel { array_t phis_2; array_t weights_2; array_t tags_2; + array_t pldis_2; npart_t offset1, offset2; + npart_t domain_idx, cntr1, cntr2; + bool use_tracking_1, use_tracking_2; const M metric; const array_t xi_min, xi_max; - const ED energy_dist; + const ED1 energy_dist_1; + const ED2 energy_dist_2; const real_t inv_V0; random_number_pool_t random_pool; @@ -57,12 +113,15 @@ namespace kernel { spidx_t spidx2, Particles& species1, Particles& species2, + npart_t inject_npart, + npart_t domain_idx, npart_t offset1, npart_t offset2, const M& metric, const array_t& xi_min, const array_t& xi_max, - const ED& energy_dist, + const ED1& energy_dist_1, + const ED2& energy_dist_2, real_t inv_V0, random_number_pool_t& random_pool) : spidx1 { spidx1 } @@ -79,6 +138,7 @@ namespace kernel { , phis_1 { species1.phi } , weights_1 { species1.weight } , tags_1 { species1.tag } + , pldis_1 { species1.pld_i } , i1s_2 { species2.i1 } , i2s_2 { species2.i2 } , i3s_2 { species2.i3 } @@ -91,26 +151,80 @@ namespace kernel { , phis_2 { species2.phi } , weights_2 { species2.weight } , tags_2 { species2.tag } + , pldis_2 { species2.pld_i } , offset1 { offset1 } , offset2 { offset2 } + , use_tracking_1 { species1.use_tracking() } + , use_tracking_2 { species2.use_tracking() } + , domain_idx { domain_idx } + , cntr1 { species1.counter() } + , cntr2 { species2.counter() } , metric { metric } , xi_min { xi_min } , xi_max { xi_max } - , energy_dist { energy_dist } + , energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } , inv_V0 { inv_V0 } - , random_pool { random_pool } {} + , random_pool { random_pool } { + if (use_tracking_1) { + printf("using tracking for species #1\n"); + species1.set_counter(cntr1 + inject_npart); +#if !defined(MPI_ENABLED) + raise::ErrorIf(species1.pld_i.extent(1) < 1, + "Particle tracking is enabled but the " + "particle integer payload size is less " + "than 1", + HERE); +#else + raise::ErrorIf(species1.pld_i.extent(1) < 2, + "Particle tracking is enabled but the " + "particle integer payload size is less " + "than 2", + HERE); +#endif + } + if (use_tracking_2) { + printf("using tracking for species #2\n"); + species2.set_counter(cntr2 + inject_npart); +#if !defined(MPI_ENABLED) + raise::ErrorIf(species2.pld_i.extent(1) < 1, + "Particle tracking is enabled but the " + "particle integer payload size is less " + "than 1", + HERE); +#else + raise::ErrorIf(species2.pld_i.extent(1) < 2, + "Particle tracking is enabled but the " + "particle integer payload size is less " + "than 2", + HERE); +#endif + } + } Inline void operator()(index_t p) const { - coord_t x_Cd { ZERO }; - vec_t v1 { ZERO }, v2 { ZERO }; + coord_t x_Cd { ZERO }; + tuple_t xi_Cd { 0 }; + tuple_t dxi_Cd { static_cast(0) }; + vec_t v1 { ZERO }, v2 { ZERO }; { // generate a random coordinate auto rand_gen = random_pool.get_state(); - x_Cd[0] = xi_min(0) + Random(rand_gen) * (xi_max(0) - xi_min(0)); + if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or + M::Dim == Dim::_3D) { + x_Cd[0] = xi_min(0) + Random(rand_gen) * (xi_max(0) - xi_min(0)); + xi_Cd[0] = static_cast(x_Cd[0]); + dxi_Cd[0] = static_cast(x_Cd[0] - xi_Cd[0]); + } if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { x_Cd[1] = xi_min(1) + Random(rand_gen) * (xi_max(1) - xi_min(1)); + xi_Cd[1] = static_cast(x_Cd[1]); + xi_Cd[1] = static_cast(x_Cd[1]); + dxi_Cd[1] = static_cast(x_Cd[1] - xi_Cd[1]); } if constexpr (M::Dim == Dim::_3D) { x_Cd[2] = xi_min(2) + Random(rand_gen) * (xi_max(2) - xi_min(2)); + xi_Cd[2] = static_cast(x_Cd[2]); + dxi_Cd[2] = static_cast(x_Cd[2] - xi_Cd[2]); } random_pool.free_state(rand_gen); } @@ -118,241 +232,72 @@ namespace kernel { coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); if constexpr (M::CoordType == Coord::Cart) { - energy_dist(x_Ph, v1, spidx1); - energy_dist(x_Ph, v2, spidx2); + energy_dist_1(x_Ph, v1, spidx1); + energy_dist_2(x_Ph, v2, spidx2); } else if constexpr (S == SimEngine::SRPIC) { coord_t x_Cd_ { ZERO }; x_Cd_[0] = x_Cd[0]; x_Cd_[1] = x_Cd[1]; x_Cd_[2] = ZERO; // phi = 0 vec_t v_Ph { ZERO }; - energy_dist(x_Ph, v_Ph, spidx1); + energy_dist_1(x_Ph, v_Ph, spidx1); metric.template transform_xyz(x_Cd_, v_Ph, v1); - energy_dist(x_Ph, v_Ph, spidx2); + energy_dist_2(x_Ph, v_Ph, spidx2); metric.template transform_xyz(x_Cd_, v_Ph, v2); } else if constexpr (S == SimEngine::GRPIC) { vec_t v_Ph { ZERO }; - energy_dist(x_Ph, v_Ph, spidx1); + energy_dist_1(x_Ph, v_Ph, spidx1); metric.template transform(x_Cd, v_Ph, v1); - energy_dist(x_Ph, v_Ph, spidx2); + energy_dist_2(x_Ph, v_Ph, spidx2); metric.template transform(x_Cd, v_Ph, v2); } else { raise::KernelError(HERE, "Unknown simulation engine"); } } - // inject - i1s_1(p + offset1) = static_cast(x_Cd[0]); - dx1s_1(p + offset1) = static_cast( - x_Cd[0] - static_cast(i1s_1(p + offset1))); - i1s_2(p + offset2) = i1s_1(p + offset1); - dx1s_2(p + offset2) = dx1s_1(p + offset1); - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - i2s_1(p + offset1) = static_cast(x_Cd[1]); - dx2s_1(p + offset1) = static_cast( - x_Cd[1] - static_cast(i2s_1(p + offset1))); - i2s_2(p + offset2) = i2s_1(p + offset1); - dx2s_2(p + offset2) = dx2s_1(p + offset1); - if constexpr (S == SimEngine::SRPIC && M::CoordType != Coord::Cart) { - phis_1(p + offset1) = ZERO; - phis_2(p + offset2) = ZERO; - } + // clang-format off + real_t weight = ONE; + if constexpr (M::CoordType != Coord::Cart) { + const auto sqrt_det_h = metric.sqrt_det_h(x_Cd); + weight = sqrt_det_h * inv_V0; } - if constexpr (M::Dim == Dim::_3D) { - i3s_1(p + offset1) = static_cast(x_Cd[2]); - dx3s_1(p + offset1) = static_cast( - x_Cd[2] - static_cast(i3s_1(p + offset1))); - i3s_2(p + offset2) = i3s_1(p + offset1); - dx3s_2(p + offset2) = dx3s_1(p + offset1); + if (not use_tracking_1) { + InjectParticle( + p + offset1, + i1s_1, i2s_1, i3s_1, + dx1s_1, dx2s_1, dx3s_1, + ux1s_1, ux2s_1, ux3s_1, + phis_1, weights_1, tags_1, pldis_1, + xi_Cd, dxi_Cd, v1, weight, ZERO); + } else { + InjectParticle( + p + offset1, + i1s_1, i2s_1, i3s_1, + dx1s_1, dx2s_1, dx3s_1, + ux1s_1, ux2s_1, ux3s_1, + phis_1, weights_1, tags_1, pldis_1, + xi_Cd, dxi_Cd, v1, weight, ZERO, domain_idx, cntr1 + p); } - ux1s_1(p + offset1) = v1[0]; - ux2s_1(p + offset1) = v1[1]; - ux3s_1(p + offset1) = v1[2]; - ux1s_2(p + offset2) = v2[0]; - ux2s_2(p + offset2) = v2[1]; - ux3s_2(p + offset2) = v2[2]; - tags_1(p + offset1) = ParticleTag::alive; - tags_2(p + offset2) = ParticleTag::alive; - if constexpr (M::CoordType == Coord::Cart) { - weights_1(p + offset1) = ONE; - weights_2(p + offset2) = ONE; + if (not use_tracking_2) { + InjectParticle( + p + offset2, + i1s_2, i2s_2, i3s_2, + dx1s_2, dx2s_2, dx3s_2, + ux1s_2, ux2s_2, ux3s_2, + phis_2, weights_2, tags_2, pldis_2, + xi_Cd, dxi_Cd, v2, weight, ZERO); } else { - const auto sqrt_det_h = metric.sqrt_det_h(x_Cd); - weights_1(p + offset1) = sqrt_det_h * inv_V0; - weights_2(p + offset2) = sqrt_det_h * inv_V0; + InjectParticle( + p + offset2, + i1s_2, i2s_2, i3s_2, + dx1s_2, dx2s_2, dx3s_2, + ux1s_2, ux2s_2, ux3s_2, + phis_2, weights_2, tags_2, pldis_2, + xi_Cd, dxi_Cd, v2, weight, ZERO, domain_idx, cntr2 + p); } + // clang-format on } }; // struct UniformInjector_kernel - namespace experimental { - - template - struct UniformInjector_kernel { - static_assert(ED1::is_energy_dist, - "ED1 must be an energy distribution class"); - static_assert(ED2::is_energy_dist, - "ED2 must be an energy distribution class"); - static_assert(M::is_metric, "M must be a metric class"); - - const spidx_t spidx1, spidx2; - - array_t i1s_1, i2s_1, i3s_1; - array_t dx1s_1, dx2s_1, dx3s_1; - array_t ux1s_1, ux2s_1, ux3s_1; - array_t phis_1; - array_t weights_1; - array_t tags_1; - - array_t i1s_2, i2s_2, i3s_2; - array_t dx1s_2, dx2s_2, dx3s_2; - array_t ux1s_2, ux2s_2, ux3s_2; - array_t phis_2; - array_t weights_2; - array_t tags_2; - - npart_t offset1, offset2; - const M metric; - const array_t xi_min, xi_max; - const ED1 energy_dist_1; - const ED2 energy_dist_2; - const real_t inv_V0; - random_number_pool_t random_pool; - - UniformInjector_kernel(spidx_t spidx1, - spidx_t spidx2, - Particles& species1, - Particles& species2, - npart_t offset1, - npart_t offset2, - const M& metric, - const array_t& xi_min, - const array_t& xi_max, - const ED1& energy_dist_1, - const ED2& energy_dist_2, - real_t inv_V0, - random_number_pool_t& random_pool) - : spidx1 { spidx1 } - , spidx2 { spidx2 } - , i1s_1 { species1.i1 } - , i2s_1 { species1.i2 } - , i3s_1 { species1.i3 } - , dx1s_1 { species1.dx1 } - , dx2s_1 { species1.dx2 } - , dx3s_1 { species1.dx3 } - , ux1s_1 { species1.ux1 } - , ux2s_1 { species1.ux2 } - , ux3s_1 { species1.ux3 } - , phis_1 { species1.phi } - , weights_1 { species1.weight } - , tags_1 { species1.tag } - , i1s_2 { species2.i1 } - , i2s_2 { species2.i2 } - , i3s_2 { species2.i3 } - , dx1s_2 { species2.dx1 } - , dx2s_2 { species2.dx2 } - , dx3s_2 { species2.dx3 } - , ux1s_2 { species2.ux1 } - , ux2s_2 { species2.ux2 } - , ux3s_2 { species2.ux3 } - , phis_2 { species2.phi } - , weights_2 { species2.weight } - , tags_2 { species2.tag } - , offset1 { offset1 } - , offset2 { offset2 } - , metric { metric } - , xi_min { xi_min } - , xi_max { xi_max } - , energy_dist_1 { energy_dist_1 } - , energy_dist_2 { energy_dist_2 } - , inv_V0 { inv_V0 } - , random_pool { random_pool } {} - - Inline void operator()(index_t p) const { - coord_t x_Cd { ZERO }; - vec_t v1 { ZERO }, v2 { ZERO }; - { // generate a random coordinate - auto rand_gen = random_pool.get_state(); - x_Cd[0] = xi_min(0) + Random(rand_gen) * (xi_max(0) - xi_min(0)); - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - x_Cd[1] = xi_min(1) + - Random(rand_gen) * (xi_max(1) - xi_min(1)); - } - if constexpr (M::Dim == Dim::_3D) { - x_Cd[2] = xi_min(2) + - Random(rand_gen) * (xi_max(2) - xi_min(2)); - } - random_pool.free_state(rand_gen); - } - { // generate the velocity - coord_t x_Ph { ZERO }; - metric.template convert(x_Cd, x_Ph); - if constexpr (M::CoordType == Coord::Cart) { - energy_dist_1(x_Ph, v1, spidx1); - energy_dist_2(x_Ph, v2, spidx2); - } else if constexpr (S == SimEngine::SRPIC) { - coord_t x_Cd_ { ZERO }; - x_Cd_[0] = x_Cd[0]; - x_Cd_[1] = x_Cd[1]; - x_Cd_[2] = ZERO; // phi = 0 - vec_t v_Ph { ZERO }; - energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform_xyz(x_Cd_, v_Ph, v1); - energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform_xyz(x_Cd_, v_Ph, v2); - } else if constexpr (S == SimEngine::GRPIC) { - vec_t v_Ph { ZERO }; - energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform(x_Cd, v_Ph, v1); - energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform(x_Cd, v_Ph, v2); - } else { - raise::KernelError(HERE, "Unknown simulation engine"); - } - } - // inject - i1s_1(p + offset1) = static_cast(x_Cd[0]); - dx1s_1(p + offset1) = static_cast( - x_Cd[0] - static_cast(i1s_1(p + offset1))); - i1s_2(p + offset2) = i1s_1(p + offset1); - dx1s_2(p + offset2) = dx1s_1(p + offset1); - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - i2s_1(p + offset1) = static_cast(x_Cd[1]); - dx2s_1(p + offset1) = static_cast( - x_Cd[1] - static_cast(i2s_1(p + offset1))); - i2s_2(p + offset2) = i2s_1(p + offset1); - dx2s_2(p + offset2) = dx2s_1(p + offset1); - if constexpr (S == SimEngine::SRPIC && M::CoordType != Coord::Cart) { - phis_1(p + offset1) = ZERO; - phis_2(p + offset2) = ZERO; - } - } - if constexpr (M::Dim == Dim::_3D) { - i3s_1(p + offset1) = static_cast(x_Cd[2]); - dx3s_1(p + offset1) = static_cast( - x_Cd[2] - static_cast(i3s_1(p + offset1))); - i3s_2(p + offset2) = i3s_1(p + offset1); - dx3s_2(p + offset2) = dx3s_1(p + offset1); - } - ux1s_1(p + offset1) = v1[0]; - ux2s_1(p + offset1) = v1[1]; - ux3s_1(p + offset1) = v1[2]; - ux1s_2(p + offset2) = v2[0]; - ux2s_2(p + offset2) = v2[1]; - ux3s_2(p + offset2) = v2[2]; - tags_1(p + offset1) = ParticleTag::alive; - tags_2(p + offset2) = ParticleTag::alive; - if constexpr (M::CoordType == Coord::Cart) { - weights_1(p + offset1) = ONE; - weights_2(p + offset2) = ONE; - } else { - const auto sqrt_det_h = metric.sqrt_det_h(x_Cd); - weights_1(p + offset1) = sqrt_det_h * inv_V0; - weights_2(p + offset2) = sqrt_det_h * inv_V0; - } - } - }; // struct UniformInjector_kernel - - } // namespace experimental - template struct GlobalInjector_kernel { static_assert(M::is_metric, "M must be a metric class"); From fa26b3278f339cb5c5a1825095d6e2ef8142148b Mon Sep 17 00:00:00 2001 From: hayk Date: Fri, 24 Oct 2025 15:50:24 -0400 Subject: [PATCH 08/15] particle output moved to prtl class --- src/framework/containers/particles.h | 10 + src/framework/containers/particles_io.cpp | 308 ++++++++++++++++++++++ src/framework/domain/output.cpp | 89 ++----- src/global/global.h | 5 + src/output/particles.h | 39 --- src/output/writer.cpp | 31 +-- src/output/writer.h | 35 ++- 7 files changed, 368 insertions(+), 149 deletions(-) delete mode 100644 src/output/particles.h diff --git a/src/framework/containers/particles.h b/src/framework/containers/particles.h index 47758b0d9..144ca611c 100644 --- a/src/framework/containers/particles.h +++ b/src/framework/containers/particles.h @@ -273,6 +273,16 @@ namespace ntt { #endif #if defined(OUTPUT_ENABLED) + void OutputDeclare(adios2::IO&) const; + + template + void OutputWrite(adios2::IO&, + adios2::Engine&, + npart_t, + std::size_t, + std::size_t, + const M&); + void CheckpointDeclare(adios2::IO&) const; void CheckpointRead(adios2::IO&, adios2::Engine&, std::size_t, std::size_t); void CheckpointWrite(adios2::IO&, adios2::Engine&, std::size_t, std::size_t) const; diff --git a/src/framework/containers/particles_io.cpp b/src/framework/containers/particles_io.cpp index 70b7f2458..d1b571d71 100644 --- a/src/framework/containers/particles_io.cpp +++ b/src/framework/containers/particles_io.cpp @@ -5,10 +5,20 @@ #include "utils/formatting.h" #include "utils/log.h" +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/minkowski.h" +#include "metrics/qkerr_schild.h" +#include "metrics/qspherical.h" +#include "metrics/spherical.h" + #include "framework/containers/particles.h" #include "output/utils/readers.h" #include "output/utils/writers.h" +#include "kernels/prtls_to_phys.hpp" + +#include #include #if defined(MPI_ENABLED) @@ -16,6 +26,275 @@ #endif namespace ntt { + /* * * * * * * * * + * Output + * * * * * * * * */ + template + void Particles::OutputDeclare(adios2::IO& io) const { + for (auto d { 0u }; d < D; ++d) { + io.DefineVariable(fmt::format("pX%d_%d", d + 1, index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + for (auto d { 0u }; d < Dim::_3D; ++d) { + io.DefineVariable(fmt::format("pU%d_%d", d + 1, index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + io.DefineVariable(fmt::format("pW_%d", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + if (npld_r() > 0) { + for (auto pr { 0 }; pr < npld_r(); ++pr) { + io.DefineVariable(fmt::format("pPLDR%d_%d", pr, index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + } + auto num_track_plds = 0; + if (use_tracking()) { +#if !defined(MPI_ENABLED) + num_track_plds = 1; + io.DefineVariable(fmt::format("pIDX_%d", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); +#else + num_track_plds = 2; + io.DefineVariable(fmt::format("pIDX_%d", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + io.DefineVariable(fmt::format("pRNK_%d", index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); +#endif + } + if (npld_i() > num_track_plds) { + for (auto pr { num_track_plds }; pr < npld_i(); ++pr) { + io.DefineVariable( + fmt::format("pPLDI%d_%d", pr - num_track_plds, index()), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } + } + } + + template + template + void Particles::OutputWrite(adios2::IO& io, + adios2::Engine& writer, + npart_t prtl_stride, + std::size_t domains_total, + std::size_t domains_offset, + const M& metric) { + if (not is_sorted()) { + RemoveDead(); + } + const npart_t nout = npart() / prtl_stride; + array_t buff_x1, buff_x2, buff_x3; + array_t buff_ux1 { "ux1", nout }; + array_t buff_ux2 { "ux2", nout }; + array_t buff_ux3 { "ux3", nout }; + array_t buff_wei { "w", nout }; + if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { + buff_x1 = array_t { "x1", nout }; + } + if constexpr (D == Dim::_2D or D == Dim::_3D) { + buff_x2 = array_t { "x2", nout }; + } + if constexpr (D == Dim::_3D or ((D == Dim::_2D) and (C != Coord::Cart))) { + buff_x3 = array_t { "x3", nout }; + } + array_t buff_pldr; + array_t buff_pldi; + + if (npld_r() > 0) { + buff_pldr = array_t { "pldr", nout, npld_r() }; + } + if (npld_i() > 0) { + buff_pldi = array_t { "pldi", nout, npld_i() }; + } + + if (nout > 0) { + // clang-format off + Kokkos::parallel_for( + "PrtlToPhys", + nout, + kernel::PrtlToPhys_kernel(prtl_stride, + buff_x1, buff_x2, buff_x3, + buff_ux1, buff_ux2, buff_ux3, + buff_wei, + buff_pldr, buff_pldi, + i1, i2, i3, + dx1, dx2, dx3, + ux1, ux2, ux3, + phi, weight, + pld_r, pld_i, + metric)); + // clang-format on + } + npart_t nout_offset = 0; + npart_t nout_total = nout; +#if defined(MPI_ENABLED) + auto nout_total_vec = std::vector(domains_total); + MPI_Allgather(&nout, + 1, + mpi::get_type(), + nout_total_vec.data(), + 1, + mpi::get_type(), + MPI_COMM_WORLD); + nout_total = 0; + for (auto r = 0; r < domains_total; ++r) { + if (r < domains_offset) { + nout_offset += nout_total_vec[r]; + } + nout_total += nout_total_vec[r]; + } +#endif // MPI_ENABLED + out::Write1DArray(io, + writer, + fmt::format("pW_%d", index()), + buff_wei, + nout, + nout_total, + nout_offset); + out::Write1DArray(io, + writer, + fmt::format("pU1_%d", index()), + buff_ux1, + nout, + nout_total, + nout_offset); + out::Write1DArray(io, + writer, + fmt::format("pU2_%d", index()), + buff_ux2, + nout, + nout_total, + nout_offset); + out::Write1DArray(io, + writer, + fmt::format("pU3_%d", index()), + buff_ux3, + nout, + nout_total, + nout_offset); + if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { + out::Write1DArray(io, + writer, + fmt::format("pX1_%d", index()), + buff_x1, + nout, + nout_total, + nout_offset); + } + if constexpr (D == Dim::_2D or D == Dim::_3D) { + out::Write1DArray(io, + writer, + fmt::format("pX2_%d", index()), + buff_x2, + nout, + nout_total, + nout_offset); + } + if constexpr (D == Dim::_3D or ((D == Dim::_2D) and (C != Coord::Cart))) { + out::Write1DArray(io, + writer, + fmt::format("pX3_%d", index()), + buff_x3, + nout, + nout_total, + nout_offset); + } + + if (npld_r() > 0) { + for (auto pr { 0 }; pr < npld_r(); ++pr) { + auto buff_sub = Kokkos::subview(buff_pldr, Kokkos::ALL, pr); + out::Write1DSubArray( + io, + writer, + fmt::format("pPLDR%d_%d", pr, index()), + buff_sub, + nout, + nout_total, + nout_offset); + } + } + auto num_track_plds = 0; + if (use_tracking()) { +#if !defined(MPI_ENABLED) + num_track_plds = 1; + { + auto buff_sub = Kokkos::subview(buff_pldi, + Kokkos::ALL, + static_cast(pldi::spcCtr)); + out::Write1DSubArray( + io, + writer, + fmt::format("pIDX_%d", pr, index()), + buff_sub, + nout, + nout_total, + nout_offset); + } +#else + num_track_plds = 2; + { + auto buff_sub = Kokkos::subview(buff_pldi, + Kokkos::ALL, + static_cast(pldi::spcCtr)); + out::Write1DSubArray( + io, + writer, + fmt::format("pIDX_%d", index()), + buff_sub, + nout, + nout_total, + nout_offset); + } + { + auto buff_sub = Kokkos::subview(buff_pldi, + Kokkos::ALL, + static_cast(pldi::domIdx)); + out::Write1DSubArray( + io, + writer, + fmt::format("pRNK_%d", index()), + buff_sub, + nout, + nout_total, + nout_offset); + } +#endif + } + if (npld_i() > num_track_plds) { + for (auto pr { num_track_plds }; pr < npld_i(); ++pr) { + auto buff_sub = Kokkos::subview(buff_pldi, + Kokkos::ALL, + static_cast(pr)); + out::Write1DSubArray( + io, + writer, + fmt::format("pPLDI%d_%d", pr - num_track_plds, index()), + buff_sub, + nout, + nout_total, + nout_offset); + } + } + } + + /* * * * * * * * * + * Checkpoints + * * * * * * * * */ template void Particles::CheckpointDeclare(adios2::IO& io) const { @@ -477,6 +756,35 @@ namespace ntt { } } +#define PARTICLES_OUTPUT_DECLARE(D, C) \ + template void Particles::OutputDeclare(adios2::IO&) const; + + PARTICLES_OUTPUT_DECLARE(Dim::_1D, Coord::Cart) + PARTICLES_OUTPUT_DECLARE(Dim::_2D, Coord::Cart) + PARTICLES_OUTPUT_DECLARE(Dim::_3D, Coord::Cart) + PARTICLES_OUTPUT_DECLARE(Dim::_2D, Coord::Sph) + PARTICLES_OUTPUT_DECLARE(Dim::_2D, Coord::Qsph) +#undef PARTICLES_OUTPUT_DECLARE + +#define PARTICLES_OUTPUT_WRITE(S, M) \ + template void Particles::OutputWrite( \ + adios2::IO&, \ + adios2::Engine&, \ + npart_t, \ + std::size_t, \ + std::size_t, \ + const M&); + + PARTICLES_OUTPUT_WRITE(SimEngine::SRPIC, metric::Minkowski) + PARTICLES_OUTPUT_WRITE(SimEngine::SRPIC, metric::Minkowski) + PARTICLES_OUTPUT_WRITE(SimEngine::SRPIC, metric::Minkowski) + PARTICLES_OUTPUT_WRITE(SimEngine::SRPIC, metric::Spherical) + PARTICLES_OUTPUT_WRITE(SimEngine::SRPIC, metric::QSpherical) + PARTICLES_OUTPUT_WRITE(SimEngine::GRPIC, metric::KerrSchild) + PARTICLES_OUTPUT_WRITE(SimEngine::GRPIC, metric::QKerrSchild) + PARTICLES_OUTPUT_WRITE(SimEngine::GRPIC, metric::KerrSchild0) +#undef PARTICLES_OUTPUT_WRITE + #define PARTICLES_CHECKPOINTS(D, C) \ template void Particles::CheckpointDeclare(adios2::IO&) const; \ template void Particles::CheckpointRead(adios2::IO&, \ diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 6dad452e8..0a41748d9 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -21,7 +21,6 @@ #include "kernels/divergences.hpp" #include "kernels/fields_to_phys.hpp" #include "kernels/particle_moments.hpp" -#include "kernels/prtls_to_phys.hpp" #include #include @@ -92,7 +91,13 @@ namespace ntt { if constexpr (M::CoordType != Coord::Cart) { dim = Dim::_3D; } - g_writer.defineParticleOutputs(dim, species_to_write); + g_writer.clearSpeciesIndex(); + for (const auto& s : species_to_write) { + g_writer.addSpeciesIndex(s); + } + for (const auto sp : g_writer.speciesIndices()) { + local_domain->species[sp - 1].OutputDeclare(g_writer.io()); + } // spectra write all particle species std::vector spectra_species {}; @@ -711,78 +716,14 @@ namespace ntt { g_writer.beginWriting(WriteMode::Particles, current_step, current_time); const auto prtl_stride = params.template get( "output.particles.stride"); - for (const auto& prtl : g_writer.speciesWriters()) { - auto& species = local_domain->species[prtl.species() - 1]; - if (not species.is_sorted()) { - species.RemoveDead(); - } - const npart_t nout = species.npart() / prtl_stride; - array_t buff_x1, buff_x2, buff_x3; - array_t buff_ux1 { "u1", nout }; - array_t buff_ux2 { "ux2", nout }; - array_t buff_ux3 { "ux3", nout }; - array_t buff_wei { "w", nout }; - if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or - M::Dim == Dim::_3D) { - buff_x1 = array_t { "x1", nout }; - } - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - buff_x2 = array_t { "x2", nout }; - } - if constexpr (M::Dim == Dim::_3D or - ((D == Dim::_2D) and (M::CoordType != Coord::Cart))) { - buff_x3 = array_t { "x3", nout }; - } - if (nout > 0) { - // clang-format off - Kokkos::parallel_for( - "PrtlToPhys", - nout, - kernel::PrtlToPhys_kernel(prtl_stride, - buff_x1, buff_x2, buff_x3, - buff_ux1, buff_ux2, buff_ux3, - buff_wei, - species.i1, species.i2, species.i3, - species.dx1, species.dx2, species.dx3, - species.ux1, species.ux2, species.ux3, - species.phi, species.weight, - local_domain->mesh.metric)); - // clang-format on - } - npart_t offset = 0; - npart_t glob_tot = nout; -#if defined(MPI_ENABLED) - auto glob_nout = std::vector(g_ndomains); - MPI_Allgather(&nout, - 1, - mpi::get_type(), - glob_nout.data(), - 1, - mpi::get_type(), - MPI_COMM_WORLD); - glob_tot = 0; - for (auto r = 0; r < g_mpi_size; ++r) { - if (r < g_mpi_rank) { - offset += glob_nout[r]; - } - glob_tot += glob_nout[r]; - } -#endif // MPI_ENABLED - g_writer.writeParticleQuantity(buff_wei, glob_tot, offset, prtl.name("W", 0)); - g_writer.writeParticleQuantity(buff_ux1, glob_tot, offset, prtl.name("U", 1)); - g_writer.writeParticleQuantity(buff_ux2, glob_tot, offset, prtl.name("U", 2)); - g_writer.writeParticleQuantity(buff_ux3, glob_tot, offset, prtl.name("U", 3)); - if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or - M::Dim == Dim::_3D) { - g_writer.writeParticleQuantity(buff_x1, glob_tot, offset, prtl.name("X", 1)); - } - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - g_writer.writeParticleQuantity(buff_x2, glob_tot, offset, prtl.name("X", 2)); - } - if constexpr (M::Dim == Dim::_3D or - ((D == Dim::_2D) and (M::CoordType != Coord::Cart))) { - g_writer.writeParticleQuantity(buff_x3, glob_tot, offset, prtl.name("X", 3)); - } + for (const auto spec : g_writer.speciesIndices()) { + local_domain->species[spec - 1].template OutputWrite( + g_writer.io(), + g_writer.writer(), + prtl_stride, + ndomains(), + local_domain->index(), + local_domain->mesh.metric); } g_writer.endWriting(WriteMode::Particles); } // end shouldWrite("particles", step, time) diff --git a/src/global/global.h b/src/global/global.h index adffcf6e9..0586a3720 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -139,6 +139,11 @@ namespace ntt { jx3 = 2 }; + enum pldi { + spcCtr = 0, + domIdx = 1 + }; + enum ParticleTag : short { dead = 0, alive diff --git a/src/output/particles.h b/src/output/particles.h deleted file mode 100644 index 0936e66f9..000000000 --- a/src/output/particles.h +++ /dev/null @@ -1,39 +0,0 @@ -/** - * @file output/particles.h - * @brief Defines the metadata for particle output - * @implements - * - out::OutputSpecies - */ - -#ifndef OUTPUT_PARTICLES_H -#define OUTPUT_PARTICLES_H - -#include "global.h" - -#include - -namespace out { - - class OutputSpecies { - const spidx_t m_sp; - - public: - OutputSpecies(spidx_t sp) : m_sp { sp } {} - - ~OutputSpecies() = default; - - [[nodiscard]] - auto species() const -> spidx_t { - return m_sp; - } - - [[nodiscard]] - auto name(const std::string& q, unsigned short c) const -> std::string { - return "p" + q + (c == 0 ? "" : std::to_string(c)) + "_" + - std::to_string(m_sp); - } - }; - -} // namespace out - -#endif // OUTPUT_PARTICLES_H diff --git a/src/output/writer.cpp b/src/output/writer.cpp index cc9ec0eb8..42d299a67 100644 --- a/src/output/writer.cpp +++ b/src/output/writer.cpp @@ -48,9 +48,8 @@ namespace out { m_trackers.insert({ type, tools::Tracker(type, interval, interval_time) }); } - auto Writer::shouldWrite(const std::string& type, - timestep_t step, - simtime_t time) -> bool { + auto Writer::shouldWrite(const std::string& type, timestep_t step, simtime_t time) + -> bool { if (m_trackers.find(type) != m_trackers.end()) { return m_trackers.at(type).shouldWrite(step, time); } else { @@ -163,32 +162,6 @@ namespace out { } } - void Writer::defineParticleOutputs(Dimension dim, - const std::vector& specs) { - m_prtl_writers.clear(); - for (const auto& s : specs) { - m_prtl_writers.emplace_back(s); - } - for (const auto& prtl : m_prtl_writers) { - for (auto d { 0u }; d < dim; ++d) { - m_io.DefineVariable(prtl.name("X", d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - } - for (auto d { 0u }; d < Dim::_3D; ++d) { - m_io.DefineVariable(prtl.name("U", d + 1), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - } - m_io.DefineVariable(prtl.name("W", 0), - { adios2::UnknownDim }, - { adios2::UnknownDim }, - { adios2::UnknownDim }); - } - } - void Writer::defineSpectraOutputs(const std::vector& specs) { m_spectra_writers.clear(); for (const auto& s : specs) { diff --git a/src/output/writer.h b/src/output/writer.h index cc3edc733..4fe0c194f 100644 --- a/src/output/writer.h +++ b/src/output/writer.h @@ -14,7 +14,6 @@ #include "utils/tools.h" #include "output/fields.h" -#include "output/particles.h" #include "output/spectra.h" #include @@ -62,9 +61,10 @@ namespace out { std::map m_trackers; std::vector m_flds_writers; - std::vector m_prtl_writers; std::vector m_spectra_writers; + std::vector m_species_indices; + WriteModeTags m_active_mode { WriteMode::None }; public: @@ -92,7 +92,6 @@ namespace out { Coord); void defineFieldOutputs(const SimEngine&, const std::vector&); - void defineParticleOutputs(Dimension, const std::vector&); void defineSpectraOutputs(const std::vector&); void writeMesh(unsigned short, @@ -115,19 +114,41 @@ namespace out { void beginWriting(WriteModeTags, timestep_t, simtime_t); void endWriting(WriteModeTags); + void addSpeciesIndex(spidx_t idx) { + m_species_indices.push_back(idx); + } + + void clearSpeciesIndex() { + m_species_indices.clear(); + } + /* getters -------------------------------------------------------------- */ + [[nodiscard]] + auto io() -> adios2::IO& { + return m_io; + } + + [[nodiscard]] + auto writer() -> adios2::Engine& { + return m_writer; + } + + [[nodiscard]] + auto speciesIndices() const -> const std::vector& { + return m_species_indices; + } + + [[nodiscard]] auto root() const -> const path_t& { return m_root; } + [[nodiscard]] auto fieldWriters() const -> const std::vector& { return m_flds_writers; } - auto speciesWriters() const -> const std::vector& { - return m_prtl_writers; - } - + [[nodiscard]] auto spectraWriters() const -> const std::vector& { return m_spectra_writers; } From fce7983847f414f1f32c562f728eb8b1455fff5e Mon Sep 17 00:00:00 2001 From: hayk Date: Fri, 24 Oct 2025 15:50:37 -0400 Subject: [PATCH 09/15] writer for subview --- src/output/utils/writers.h | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/output/utils/writers.h b/src/output/utils/writers.h index 58fac2bf6..1da24e9ee 100644 --- a/src/output/utils/writers.h +++ b/src/output/utils/writers.h @@ -41,6 +41,33 @@ namespace out { std::size_t, std::size_t); + // template + // void Write1DSubArray(adios2::IO&, + // adios2::Engine&, + // const std::string&, + // const subarray1d_t&, + // std::size_t, + // std::size_t, + // std::size_t); + + template + void Write1DSubArray(adios2::IO& io, + adios2::Engine& writer, + const std::string& name, + const S& data, + std::size_t local_size, + std::size_t global_size, + std::size_t local_offset) { + const auto slice = range_tuple_t(0, local_size); + auto var = io.InquireVariable(name); + var.SetShape({ global_size }); + var.SetSelection(adios2::Box({ local_offset }, { local_size })); + + auto data_h = Kokkos::create_mirror_view(data); + Kokkos::deep_copy(data_h, data); + writer.Put(var, data_h.data(), adios2::Mode::Sync); + } + template void Write2DArray(adios2::IO&, adios2::Engine&, From 3a35b75320ef21903b9165e025227a8f5bf8a4bb Mon Sep 17 00:00:00 2001 From: hayk Date: Fri, 24 Oct 2025 15:50:56 -0400 Subject: [PATCH 10/15] pld passed to prtl2phys kernel --- src/kernels/prtls_to_phys.hpp | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/kernels/prtls_to_phys.hpp b/src/kernels/prtls_to_phys.hpp index 4dd7d88b0..dfe039e86 100644 --- a/src/kernels/prtls_to_phys.hpp +++ b/src/kernels/prtls_to_phys.hpp @@ -39,11 +39,15 @@ namespace kernel { array_t buff_ux2; array_t buff_ux3; array_t buff_wei; + array_t buff_pldr; + array_t buff_pldi; const array_t i1, i2, i3; const array_t dx1, dx2, dx3; const array_t ux1, ux2, ux3; const array_t phi; const array_t weight; + const array_t pld_r; + const array_t pld_i; const M metric; public: @@ -55,6 +59,8 @@ namespace kernel { array_t& buff_ux2, array_t& buff_ux3, array_t& buff_wei, + array_t& buff_pldr, + array_t& buff_pldi, const array_t& i1, const array_t& i2, const array_t& i3, @@ -66,6 +72,8 @@ namespace kernel { const array_t& ux3, const array_t& phi, const array_t& weight, + const array_t& pld_r, + const array_t& pld_i, const M& metric) : stride { stride } , buff_x1 { buff_x1 } @@ -75,6 +83,8 @@ namespace kernel { , buff_ux2 { buff_ux2 } , buff_ux3 { buff_ux3 } , buff_wei { buff_wei } + , buff_pldr { buff_pldr } + , buff_pldi { buff_pldi } , i1 { i1 } , i2 { i2 } , i3 { i3 } @@ -86,6 +96,8 @@ namespace kernel { , ux3 { ux3 } , phi { phi } , weight { weight } + , pld_r { pld_r } + , pld_i { pld_i } , metric { metric } { if constexpr ((D == Dim::_1D) || (D == Dim::_2D) || (D == Dim::_3D)) { raise::ErrorIf(buff_x1.extent(0) == 0, "Invalid buffer size", HERE); @@ -106,6 +118,7 @@ namespace kernel { bufferX(p); bufferU(p); buff_wei(p) = weight(p * stride); + bufferPlds(p); } Inline void bufferX(index_t& p) const { @@ -199,6 +212,19 @@ namespace kernel { buff_ux2(p) = u_Phys[1]; buff_ux3(p) = u_Phys[2]; } + + Inline void bufferPlds(index_t& p) const { + if (buff_pldr.extent(0) > 0) { + for (auto pr { 0u }; pr < buff_pldr.extent(1); ++pr) { + buff_pldr(p, pr) = pld_r(p * stride, pr); + } + } + if (buff_pldi.extent(0) > 0) { + for (auto pi { 0u }; pi < buff_pldi.extent(1); ++pi) { + buff_pldi(p, pi) = pld_i(p * stride, pi); + } + } + } }; } // namespace kernel From bb41984bcbd2a0fcc5d5787beb06b5ba3b7d996c Mon Sep 17 00:00:00 2001 From: hayk Date: Fri, 24 Oct 2025 15:51:10 -0400 Subject: [PATCH 11/15] minor formatting --- src/global/arch/kokkos_aliases.h | 4 ++-- src/kernels/particle_pusher_sr.hpp | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/global/arch/kokkos_aliases.h b/src/global/arch/kokkos_aliases.h index adb0b6451..712fc6eff 100644 --- a/src/global/arch/kokkos_aliases.h +++ b/src/global/arch/kokkos_aliases.h @@ -234,8 +234,8 @@ auto CreateParticleRangePolicy(npart_t, npart_t) -> range_t; * @returns Kokkos::RangePolicy or Kokkos::MDRangePolicy in the accelerator execution space. */ template -auto CreateRangePolicy(const tuple_t&, - const tuple_t&) -> range_t; +auto CreateRangePolicy(const tuple_t&, const tuple_t&) + -> range_t; /** * @brief Function template for generating ND Kokkos range policy on the host. diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 91bc6a760..6bd4e1714 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -30,7 +30,9 @@ /* Local macros */ /* -------------------------------------------------------------------------- */ #define from_Xi_to_i(XI, I) \ - { I = static_cast((XI + 1)) - 1; } + { \ + I = static_cast((XI + 1)) - 1; \ + } #define from_Xi_to_i_di(XI, I, DI) \ { \ From 713fd6e585326edc0b8e1eb20d7ffbc3e4b142bf Mon Sep 17 00:00:00 2001 From: haykh Date: Mon, 27 Oct 2025 10:52:54 -0400 Subject: [PATCH 12/15] prtl tracking in output --- src/framework/containers/particles_io.cpp | 119 +++++++++++++++------- src/kernels/prtls_to_phys.hpp | 106 +++++++++---------- 2 files changed, 134 insertions(+), 91 deletions(-) diff --git a/src/framework/containers/particles_io.cpp b/src/framework/containers/particles_io.cpp index d1b571d71..9131e596a 100644 --- a/src/framework/containers/particles_io.cpp +++ b/src/framework/containers/particles_io.cpp @@ -97,7 +97,56 @@ namespace ntt { if (not is_sorted()) { RemoveDead(); } - const npart_t nout = npart() / prtl_stride; + npart_t nout; + array_t out_indices; + if (!use_tracking()) { + nout = npart() / prtl_stride; + } else { + nout = 0u; + Kokkos::parallel_reduce( + "CountOutputParticles", + npart(), + Lambda(index_t p, npart_t & l_nout) { + if ((tag(p) == ParticleTag::alive) and + (pld_i(p, pldi::spcCtr) % prtl_stride == 0)) { + l_nout += 1; + } + }, + nout); + out_indices = array_t { "out_indices", nout }; + array_t out_counter { "out_counter" }; + Kokkos::parallel_for( + "RecordOutputIndices", + npart(), + Lambda(index_t p) { + if ((tag(p) == ParticleTag::alive) and + (pld_i(p, pldi::spcCtr) % prtl_stride == 0)) { + const auto p_out = Kokkos::atomic_fetch_add(&out_counter(), 1); + out_indices(p_out) = p; + } + }); + } + + npart_t nout_offset = 0; + npart_t nout_total = nout; +#if defined(MPI_ENABLED) + auto nout_total_vec = std::vector(domains_total); + MPI_Allgather(&nout, + 1, + mpi::get_type(), + nout_total_vec.data(), + 1, + mpi::get_type(), + MPI_COMM_WORLD); + nout_total = 0; + for (auto r = 0; r < domains_total; ++r) { + if (r < domains_offset) { + nout_offset += nout_total_vec[r]; + } + nout_total += nout_total_vec[r]; + } +#endif // MPI_ENABLED + array_t buff_x1, buff_x2, buff_x3; array_t buff_ux1 { "ux1", nout }; array_t buff_ux2 { "ux2", nout }; @@ -123,42 +172,42 @@ namespace ntt { } if (nout > 0) { - // clang-format off - Kokkos::parallel_for( - "PrtlToPhys", - nout, - kernel::PrtlToPhys_kernel(prtl_stride, - buff_x1, buff_x2, buff_x3, - buff_ux1, buff_ux2, buff_ux3, - buff_wei, - buff_pldr, buff_pldi, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, - pld_r, pld_i, - metric)); - // clang-format on - } - npart_t nout_offset = 0; - npart_t nout_total = nout; -#if defined(MPI_ENABLED) - auto nout_total_vec = std::vector(domains_total); - MPI_Allgather(&nout, - 1, - mpi::get_type(), - nout_total_vec.data(), - 1, - mpi::get_type(), - MPI_COMM_WORLD); - nout_total = 0; - for (auto r = 0; r < domains_total; ++r) { - if (r < domains_offset) { - nout_offset += nout_total_vec[r]; + if (!use_tracking()) { + // clang-format off + Kokkos::parallel_for( + "PrtlToPhys", + nout, + kernel::PrtlToPhys_kernel(prtl_stride, out_indices, + buff_x1, buff_x2, buff_x3, + buff_ux1, buff_ux2, buff_ux3, + buff_wei, + buff_pldr, buff_pldi, + i1, i2, i3, + dx1, dx2, dx3, + ux1, ux2, ux3, + phi, weight, + pld_r, pld_i, + metric)); + // clang-format on + } else { + // clang-format off + Kokkos::parallel_for( + "PrtlToPhys", + nout, + kernel::PrtlToPhys_kernel(prtl_stride, out_indices, + buff_x1, buff_x2, buff_x3, + buff_ux1, buff_ux2, buff_ux3, + buff_wei, + buff_pldr, buff_pldi, + i1, i2, i3, + dx1, dx2, dx3, + ux1, ux2, ux3, + phi, weight, + pld_r, pld_i, + metric)); + // clang-format on } - nout_total += nout_total_vec[r]; } -#endif // MPI_ENABLED out::Write1DArray(io, writer, fmt::format("pW_%d", index()), diff --git a/src/kernels/prtls_to_phys.hpp b/src/kernels/prtls_to_phys.hpp index dfe039e86..fbafbe00e 100644 --- a/src/kernels/prtls_to_phys.hpp +++ b/src/kernels/prtls_to_phys.hpp @@ -25,13 +25,14 @@ namespace kernel { using namespace ntt; - template + template class PrtlToPhys_kernel { static_assert(M::is_metric, "M must be a metric class"); static constexpr Dimension D = M::Dim; protected: const npart_t stride; + array_t out_indices; array_t buff_x1; array_t buff_x2; array_t buff_x3; @@ -52,6 +53,7 @@ namespace kernel { public: PrtlToPhys_kernel(npart_t stride, + array_t out_indices, array_t& buff_x1, array_t& buff_x2, array_t& buff_x3, @@ -76,6 +78,7 @@ namespace kernel { const array_t& pld_i, const M& metric) : stride { stride } + , out_indices { out_indices } , buff_x1 { buff_x1 } , buff_x2 { buff_x2 } , buff_x3 { buff_x3 } @@ -115,41 +118,44 @@ namespace kernel { } Inline void operator()(index_t p) const { - bufferX(p); - bufferU(p); - buff_wei(p) = weight(p * stride); - bufferPlds(p); + if constexpr (!T) { // no tracking enabled + bufferX(p * stride, p); + bufferU(p * stride, p); + buff_wei(p) = weight(p * stride); + bufferPlds(p * stride, p); + } else { + bufferX(out_indices(p), p); + bufferU(out_indices(p), p); + buff_wei(p) = weight(out_indices(p)); + bufferPlds(out_indices(p), p); + } } - Inline void bufferX(index_t& p) const { + Inline void bufferX(index_t& p_from, index_t& p_to) const { if constexpr ((D == Dim::_1D) || (D == Dim::_2D) || (D == Dim::_3D)) { - buff_x1(p) = metric.template convert<1, Crd::Cd, Crd::Ph>( - static_cast(i1(p * stride)) + - static_cast(dx1(p * stride))); + buff_x1(p_to) = metric.template convert<1, Crd::Cd, Crd::Ph>( + static_cast(i1(p_from)) + static_cast(dx1(p_from))); } if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) { - buff_x2(p) = metric.template convert<2, Crd::Cd, Crd::Ph>( - static_cast(i2(p * stride)) + - static_cast(dx2(p * stride))); + buff_x2(p_to) = metric.template convert<2, Crd::Cd, Crd::Ph>( + static_cast(i2(p_from)) + static_cast(dx2(p_from))); } if constexpr ((D == Dim::_2D) && (M::CoordType != Coord::Cart)) { - buff_x3(p) = phi(p * stride); + buff_x3(p_to) = phi(p_from); } if constexpr (D == Dim::_3D) { - buff_x3(p) = metric.template convert<3, Crd::Cd, Crd::Ph>( - static_cast(i3(p * stride)) + - static_cast(dx3(p * stride))); + buff_x3(p_to) = metric.template convert<3, Crd::Cd, Crd::Ph>( + static_cast(i3(p_from)) + static_cast(dx3(p_from))); } } - Inline void bufferU(index_t& p) const { + Inline void bufferU(index_t& p_from, index_t& p_to) const { vec_t u_Phys { ZERO }; if constexpr (D == Dim::_1D) { if constexpr (M::CoordType == Coord::Cart) { metric.template transform_xyz( - { static_cast(i1(p * stride)) + - static_cast(dx1(p * stride)) }, - { ux1(p * stride), ux2(p * stride), ux3(p * stride) }, + { static_cast(i1(p_from)) + static_cast(dx1(p_from)) }, + { ux1(p_from), ux2(p_from), ux3(p_from) }, u_Phys); } else { raise::KernelError(HERE, "Unsupported coordinate system in 1D"); @@ -157,28 +163,22 @@ namespace kernel { } else if constexpr (D == Dim::_2D) { if constexpr (M::CoordType == Coord::Cart) { metric.template transform_xyz( - { static_cast(i1(p * stride)) + - static_cast(dx1(p * stride)), - static_cast(i2(p * stride)) + - static_cast(dx2(p * stride)) }, - { ux1(p * stride), ux2(p * stride), ux3(p * stride) }, + { static_cast(i1(p_from)) + static_cast(dx1(p_from)), + static_cast(i2(p_from)) + static_cast(dx2(p_from)) }, + { ux1(p_from), ux2(p_from), ux3(p_from) }, u_Phys); } else if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz( - { static_cast(i1(p * stride)) + - static_cast(dx1(p * stride)), - static_cast(i2(p * stride)) + - static_cast(dx2(p * stride)), - phi(p * stride) }, - { ux1(p * stride), ux2(p * stride), ux3(p * stride) }, + { static_cast(i1(p_from)) + static_cast(dx1(p_from)), + static_cast(i2(p_from)) + static_cast(dx2(p_from)), + phi(p_from) }, + { ux1(p_from), ux2(p_from), ux3(p_from) }, u_Phys); } else if constexpr (S == SimEngine::GRPIC) { metric.template transform( - { static_cast(i1(p * stride)) + - static_cast(dx1(p * stride)), - static_cast(i2(p * stride)) + - static_cast(dx2(p * stride)) }, - { ux1(p * stride), ux2(p * stride), ux3(p * stride) }, + { static_cast(i1(p_from)) + static_cast(dx1(p_from)), + static_cast(i2(p_from)) + static_cast(dx2(p_from)) }, + { ux1(p_from), ux2(p_from), ux3(p_from) }, u_Phys); } else { raise::KernelError(HERE, "Unrecognized simulation engine"); @@ -186,42 +186,36 @@ namespace kernel { } else if constexpr (D == Dim::_3D) { if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz( - { static_cast(i1(p * stride)) + - static_cast(dx1(p * stride)), - static_cast(i2(p * stride)) + - static_cast(dx2(p * stride)), - static_cast(i3(p * stride)) + - static_cast(dx3(p * stride)) }, - { ux1(p * stride), ux2(p * stride), ux3(p * stride) }, + { static_cast(i1(p_from)) + static_cast(dx1(p_from)), + static_cast(i2(p_from)) + static_cast(dx2(p_from)), + static_cast(i3(p_from)) + static_cast(dx3(p_from)) }, + { ux1(p_from), ux2(p_from), ux3(p_from) }, u_Phys); } else if constexpr (S == SimEngine::GRPIC) { metric.template transform( - { static_cast(i1(p * stride)) + - static_cast(dx1(p * stride)), - static_cast(i2(p * stride)) + - static_cast(dx2(p * stride)), - static_cast(i3(p * stride)) + - static_cast(dx3(p * stride)) }, - { ux1(p * stride), ux2(p * stride), ux3(p * stride) }, + { static_cast(i1(p_from)) + static_cast(dx1(p_from)), + static_cast(i2(p_from)) + static_cast(dx2(p_from)), + static_cast(i3(p_from)) + static_cast(dx3(p_from)) }, + { ux1(p_from), ux2(p_from), ux3(p_from) }, u_Phys); } else { raise::KernelError(HERE, "Unrecognized simulation engine"); } } - buff_ux1(p) = u_Phys[0]; - buff_ux2(p) = u_Phys[1]; - buff_ux3(p) = u_Phys[2]; + buff_ux1(p_to) = u_Phys[0]; + buff_ux2(p_to) = u_Phys[1]; + buff_ux3(p_to) = u_Phys[2]; } - Inline void bufferPlds(index_t& p) const { + Inline void bufferPlds(index_t& p_from, index_t& p_to) const { if (buff_pldr.extent(0) > 0) { for (auto pr { 0u }; pr < buff_pldr.extent(1); ++pr) { - buff_pldr(p, pr) = pld_r(p * stride, pr); + buff_pldr(p_to, pr) = pld_r(p_from, pr); } } if (buff_pldi.extent(0) > 0) { for (auto pi { 0u }; pi < buff_pldi.extent(1); ++pi) { - buff_pldi(p, pi) = pld_i(p * stride, pi); + buff_pldi(p_to, pi) = pld_i(p_from, pi); } } } From dbc039f3f1c439bcca66de2de7bb6623bbfa4ac0 Mon Sep 17 00:00:00 2001 From: haykh Date: Fri, 7 Nov 2025 10:09:29 -0800 Subject: [PATCH 13/15] capture this explicitly --- src/framework/containers/particles_io.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/framework/containers/particles_io.cpp b/src/framework/containers/particles_io.cpp index 9131e596a..332c6e9d9 100644 --- a/src/framework/containers/particles_io.cpp +++ b/src/framework/containers/particles_io.cpp @@ -102,13 +102,15 @@ namespace ntt { if (!use_tracking()) { nout = npart() / prtl_stride; } else { - nout = 0u; + nout = 0u; + const auto tag_d = this->tag; + const auto pld_i_d = this->pld_i; Kokkos::parallel_reduce( "CountOutputParticles", npart(), Lambda(index_t p, npart_t & l_nout) { - if ((tag(p) == ParticleTag::alive) and - (pld_i(p, pldi::spcCtr) % prtl_stride == 0)) { + if ((tag_d(p) == ParticleTag::alive) and + (pld_i_d(p, pldi::spcCtr) % prtl_stride == 0)) { l_nout += 1; } }, @@ -119,8 +121,8 @@ namespace ntt { "RecordOutputIndices", npart(), Lambda(index_t p) { - if ((tag(p) == ParticleTag::alive) and - (pld_i(p, pldi::spcCtr) % prtl_stride == 0)) { + if ((tag_d(p) == ParticleTag::alive) and + (pld_i_d(p, pldi::spcCtr) % prtl_stride == 0)) { const auto p_out = Kokkos::atomic_fetch_add(&out_counter(), 1); out_indices(p_out) = p; } @@ -288,7 +290,7 @@ namespace ntt { out::Write1DSubArray( io, writer, - fmt::format("pIDX_%d", pr, index()), + fmt::format("pIDX_%d", index()), buff_sub, nout, nout_total, From 2e367b0f99771f3d85848dbecad229aef60b5dfe Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 8 Nov 2025 09:19:34 -0800 Subject: [PATCH 14/15] tests adjusted --- cmake/tests.cmake | 7 -- src/kernels/tests/prtls_to_phys.cpp | 103 +++++++++++++++++----------- src/output/tests/writer-mpi.cpp | 31 +++++---- 3 files changed, 80 insertions(+), 61 deletions(-) diff --git a/cmake/tests.cmake b/cmake/tests.cmake index 189cc2cc4..0eb043f70 100644 --- a/cmake/tests.cmake +++ b/cmake/tests.cmake @@ -9,9 +9,6 @@ add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) add_subdirectory(${SRC_DIR}/output ${CMAKE_CURRENT_BINARY_DIR}/output) -if(${output}) - add_subdirectory(${SRC_DIR}/checkpoint ${CMAKE_CURRENT_BINARY_DIR}/checkpoint) -endif() set(TEST_DIRECTORIES "") @@ -27,10 +24,6 @@ endif() list(APPEND TEST_DIRECTORIES output) -if(${output}) - list(APPEND TEST_DIRECTORIES checkpoint) -endif() - foreach(test_dir IN LISTS TEST_DIRECTORIES) add_subdirectory(${SRC_DIR}/${test_dir}/tests ${CMAKE_CURRENT_BINARY_DIR}/${test_dir}/tests) diff --git a/src/kernels/tests/prtls_to_phys.cpp b/src/kernels/tests/prtls_to_phys.cpp index 962c21b5c..8f9d1760a 100644 --- a/src/kernels/tests/prtls_to_phys.cpp +++ b/src/kernels/tests/prtls_to_phys.cpp @@ -19,9 +19,7 @@ #include #include -#include #include -#include #include using namespace ntt; @@ -39,13 +37,15 @@ struct Checker { const array_t& ux2, const array_t& ux3, const array_t& weight, + const array_t& pld_i, const array_t& buff_x1, const array_t& buff_x2, const array_t& buff_x3, const array_t& buff_ux1, const array_t& buff_ux2, const array_t& buff_ux3, - const array_t& buff_wei) + const array_t& buff_wei, + const array_t& buff_pld_i) : metric { metric } , stride { stride } , i1 { i1 } @@ -57,13 +57,15 @@ struct Checker { , ux2 { ux2 } , ux3 { ux3 } , weight { weight } + , pld_i { pld_i } , buff_x1 { buff_x1 } , buff_x2 { buff_x2 } , buff_x3 { buff_x3 } , buff_ux1 { buff_ux1 } , buff_ux2 { buff_ux2 } , buff_ux3 { buff_ux3 } - , buff_wei { buff_wei } {} + , buff_wei { buff_wei } + , buff_pld_i { buff_pld_i } {} Inline void operator()(index_t p) const { std::size_t pold = p * stride; @@ -97,6 +99,9 @@ struct Checker { if (not cmp::AlmostEqual(weight(pold), buff_wei(p))) { raise::KernelError(HERE, "weight != buff_wei"); } + if (pld_i(pold, pldi::spcCtr) != buff_pld_i(p, pldi::spcCtr)) { + raise::KernelError(HERE, "weight != buff_wei"); + } } private: @@ -111,13 +116,15 @@ struct Checker { const array_t ux2; const array_t ux3; const array_t weight; - array_t buff_x1; - array_t buff_x2; - array_t buff_x3; - array_t buff_ux1; - array_t buff_ux2; - array_t buff_ux3; - array_t buff_wei; + const array_t pld_i; + const array_t buff_x1; + const array_t buff_x2; + const array_t buff_x3; + const array_t buff_ux1; + const array_t buff_ux2; + const array_t buff_ux3; + const array_t buff_wei; + const array_t buff_pld_i; }; template @@ -148,10 +155,14 @@ void testPrtl2PhysSR(const std::vector& res, array_t ux2 { "ux2", nprtl }; array_t ux3 { "ux3", nprtl }; array_t weight { "weight", nprtl }; + array_t pldr; + array_t pld_i { "pld_i", nprtl, 1 }; array_t i3; array_t dx3; + const std::size_t stride = 2; + array_t out_indices { "out_indices", nprtl / stride }; Kokkos::parallel_for( "Init", nprtl, @@ -166,40 +177,50 @@ void testPrtl2PhysSR(const std::vector& res, ux2(p) = ((real_t)(p) - (real_t)(nprtl) / 4) / (real_t)(9 * nprtl); ux3(p) = ((real_t)(p) - (real_t)(nprtl) / 2) / (real_t)(5 * nprtl); weight(p) = (real_t)(25) + (real_t)(p) / (real_t)(nprtl); + pld_i(p, pldi::spcCtr) = p; + if (p % stride == 0) { + out_indices(p / stride) = p; + } }); - const std::size_t stride = 2; - array_t buff_x1 { "buff_x1", nprtl / stride }; - array_t buff_x2 { "buff_x2", nprtl / stride }; - array_t buff_x3 { "buff_x3", nprtl / stride }; - array_t buff_ux1 { "buff_ux1", nprtl / stride }; - array_t buff_ux2 { "buff_ux2", nprtl / stride }; - array_t buff_ux3 { "buff_ux3", nprtl / stride }; - array_t buff_wei { "buff_wei", nprtl / stride }; + array_t buff_x1 { "buff_x1", nprtl / stride }; + array_t buff_x2 { "buff_x2", nprtl / stride }; + array_t buff_x3 { "buff_x3", nprtl / stride }; + array_t buff_ux1 { "buff_ux1", nprtl / stride }; + array_t buff_ux2 { "buff_ux2", nprtl / stride }; + array_t buff_ux3 { "buff_ux3", nprtl / stride }; + array_t buff_wei { "buff_wei", nprtl / stride }; + array_t buff_pldr; + array_t buff_pld_i { "pld_i", nprtl / stride, 1 }; Kokkos::parallel_for( "Init", Kokkos::RangePolicy(0, nprtl / stride), - kernel::PrtlToPhys_kernel(stride, - buff_x1, - buff_x2, - buff_x3, - buff_ux1, - buff_ux2, - buff_ux3, - buff_wei, - i1, - i2, - i3, - dx1, - dx2, - dx3, - ux1, - ux2, - ux3, - phi, - weight, - metric)); + kernel::PrtlToPhys_kernel(stride, + out_indices, + buff_x1, + buff_x2, + buff_x3, + buff_ux1, + buff_ux2, + buff_ux3, + buff_wei, + buff_pldr, + buff_pld_i, + i1, + i2, + i3, + dx1, + dx2, + dx3, + ux1, + ux2, + ux3, + phi, + weight, + pldr, + pld_i, + metric)); Kokkos::parallel_for("Check", nprtl / stride, Checker(metric, @@ -213,13 +234,15 @@ void testPrtl2PhysSR(const std::vector& res, ux2, ux3, weight, + pld_i, buff_x1, buff_x2, buff_x3, buff_ux1, buff_ux2, buff_ux3, - buff_wei)); + buff_wei, + buff_pld_i)); } auto main(int argc, char* argv[]) -> int { diff --git a/src/output/tests/writer-mpi.cpp b/src/output/tests/writer-mpi.cpp index bc95bbc81..dad981a40 100644 --- a/src/output/tests/writer-mpi.cpp +++ b/src/output/tests/writer-mpi.cpp @@ -17,8 +17,8 @@ void cleanup() { namespace fs = std::filesystem; - fs::path tempfile_path { "test.h5" }; - fs::remove(tempfile_path); + fs::path tempfile_path { "test.bp" }; + fs::remove_all(tempfile_path); } #define CEILDIV(a, b) \ @@ -61,7 +61,7 @@ auto main(int argc, char* argv[]) -> int { { // write auto writer = out::Writer(); - writer.init(&adios, "hdf5", "test", false); + writer.init(&adios, "BPFile", "test", false); writer.defineMeshLayout({ static_cast(mpi_size) * nx1 }, { static_cast(mpi_rank) * nx1 }, { nx1 }, @@ -91,16 +91,16 @@ auto main(int argc, char* argv[]) -> int { { // read adios2::IO io = adios.DeclareIO("read-test"); - io.SetEngine("HDF5"); - adios2::Engine reader = io.Open("test.h5", adios2::Mode::Read); - raise::ErrorIf(io.InquireAttribute("NGhosts").Data()[0] != 0, - "NGhosts is not correct", - HERE); - raise::ErrorIf(io.InquireAttribute("Dimension").Data()[0] != 1, - "Dimension is not correct", - HERE); - for (std::size_t step = 0; reader.BeginStep() == adios2::StepStatus::OK; - ++step) { + io.SetEngine("BPFile"); + adios2::Engine reader = io.Open("test.bp", adios2::Mode::Read); + for (auto step = 0u; reader.BeginStep() == adios2::StepStatus::OK; ++step) { + raise::ErrorIf(io.InquireAttribute("NGhosts").Data()[0] != 0, + "NGhosts is not correct", + HERE); + raise::ErrorIf(io.InquireAttribute("Dimension").Data()[0] != 1, + "Dimension is not correct", + HERE); + timestep_t step_read; simtime_t time_read; @@ -173,6 +173,7 @@ auto main(int argc, char* argv[]) -> int { } ++cntr; } + reader.EndStep(); } reader.Close(); } @@ -186,7 +187,9 @@ auto main(int argc, char* argv[]) -> int { Kokkos::finalize(); return 1; } - cleanup(); + CallOnce([]() { + cleanup(); + }); MPI_Finalize(); Kokkos::finalize(); return 0; From ca093bcc9a8e71397cf786d6367e2605d3d668f8 Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 8 Nov 2025 09:19:43 -0800 Subject: [PATCH 15/15] bump dev/nix versions --- dev/nix/adios2.nix | 18 ++++++++++++++++-- dev/nix/kokkos.nix | 15 ++------------- dev/nix/shell.nix | 43 +++++++++++++++++++++---------------------- 3 files changed, 39 insertions(+), 37 deletions(-) diff --git a/dev/nix/adios2.nix b/dev/nix/adios2.nix index fb574c302..a3890b788 100644 --- a/dev/nix/adios2.nix +++ b/dev/nix/adios2.nix @@ -43,9 +43,23 @@ stdenv.mkDerivation { ] ++ ( if hdf5 then - (if mpi then [ pkgs.hdf5-mpi ] else [ pkgs.hdf5-cpp ]) + ( + if mpi then + [ + pkgs.hdf5-mpi + ] + else + [ pkgs.hdf5-cpp ] + ) else - (if mpi then [ pkgs.mpi ] else [ ]) + ( + if mpi then + [ + pkgs.mpi + ] + else + [ ] + ) ); configurePhase = '' diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix index 2f6ee6b99..7d86e665b 100644 --- a/dev/nix/kokkos.nix +++ b/dev/nix/kokkos.nix @@ -7,7 +7,7 @@ let name = "kokkos"; - pversion = "4.6.01"; + pversion = "4.7.01"; compilerPkgs = { "HIP" = with pkgs.rocmPackages; [ llvm.rocm-merged-llvm @@ -56,7 +56,7 @@ pkgs.stdenv.mkDerivation rec { src = pkgs.fetchgit { url = "https://github.com/kokkos/kokkos/"; rev = "${pversion}"; - sha256 = "sha256-+yszUbdHqhIkJZiGLZ9Ln4DYUosuJWKhO8FkbrY0/tY="; + sha256 = "sha256-MgphOsKE8umgYxVQZzex+elgvDDC09JaMCoU5YXaLco="; }; nativeBuildInputs = with pkgs; [ @@ -92,15 +92,4 @@ pkgs.stdenv.mkDerivation rec { installPhase = '' cmake --install build ''; - - # cmakeFlags = [ - # "-D CMAKE_CXX_STANDARD=17" - # "-D CMAKE_CXX_EXTENSIONS=OFF" - # "-D CMAKE_POSITION_INDEPENDENT_CODE=TRUE" - # "-D Kokkos_ARCH_${getArch { }}=ON" - # (if gpu != "none" then "-D Kokkos_ENABLE_${gpu}=ON" else "") - # "-D CMAKE_BUILD_TYPE=Release" - # ] ++ (cmakeExtraFlags.${gpu} src); - - # enableParallelBuilding = true; } diff --git a/dev/nix/shell.nix b/dev/nix/shell.nix index 33ae57095..0d4cc9119 100644 --- a/dev/nix/shell.nix +++ b/dev/nix/shell.nix @@ -5,7 +5,7 @@ }, gpu ? "NONE", arch ? "NATIVE", - hdf5 ? true, + hdf5 ? false, mpi ? false, }: @@ -62,27 +62,26 @@ pkgs.mkShell { pkgs.zlib ]); - shellHook = - '' - BLUE='\033[0;34m' - NC='\033[0m' + shellHook = '' + BLUE='\033[0;34m' + NC='\033[0m' - echo "following environment variables are set:" - '' - + pkgs.lib.concatStringsSep "" ( - pkgs.lib.mapAttrsToList ( - category: vars: - pkgs.lib.concatStringsSep "" ( - pkgs.lib.mapAttrsToList (name: value: '' - export ${name}=${value} - echo -e " ''\${BLUE}${name}''\${NC}=${value}" - '') vars.${gpuUpper} - ) - ) envVars - ) - + '' - echo "" - echo -e "${name} nix-shell activated" - ''; + echo "following environment variables are set:" + '' + + pkgs.lib.concatStringsSep "" ( + pkgs.lib.mapAttrsToList ( + category: vars: + pkgs.lib.concatStringsSep "" ( + pkgs.lib.mapAttrsToList (name: value: '' + export ${name}=${value} + echo -e " ''\${BLUE}${name}''\${NC}=${value}" + '') vars.${gpuUpper} + ) + ) envVars + ) + + '' + echo "" + echo -e "${name} nix-shell activated" + ''; }