diff --git a/.github/workflows/ci-macos.yml b/.github/workflows/ci-macos.yml index a1666b5a9520..45edb719c746 100644 --- a/.github/workflows/ci-macos.yml +++ b/.github/workflows/ci-macos.yml @@ -28,10 +28,14 @@ jobs: cache-dependency-path: '**/requirements.txt' - run: pip install -r requirements.txt - - name: Install dependencies + - name: Install dependencies (Homebrew) run: | brew install openmpi hdf5-mpi adios2 || true + - name: Install OpenPMD + run: | + openPMD_USE_MPI=ON python3 -m pip install openpmd-api --no-binary openpmd-api + - name: Configure run: cmake -B build -DCMAKE_BUILD_TYPE=Release diff --git a/CMakeLists.txt b/CMakeLists.txt index 590f277ee8d5..b742b21df5ee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,7 @@ option(PARTHENON_DISABLE_MPI "MPI is enabled by default if found, set this to Tr option(PARTHENON_ENABLE_HOST_COMM_BUFFERS "CUDA/HIP Only: Allocate communication buffers on host (may be slower)" OFF) option(PARTHENON_DISABLE_HDF5 "HDF5 is enabled by default if found, set this to True to disable HDF5" OFF) option(PARTHENON_DISABLE_HDF5_COMPRESSION "HDF5 compression is enabled by default, set this to True to disable compression in HDF5 output/restart files" OFF) +option(PARTHENON_ENABLE_OPENPMD "OpenPMD is enabled by default if found, set this to True to disable OpenPMD" ON) option(PARTHENON_DISABLE_SPARSE "Sparse capability is enabled by default, set this to True to compile-time disable all sparse capability" OFF) option(PARTHENON_ENABLE_ASCENT "Enable Ascent for in situ visualization and analysis" OFF) option(PARTHENON_LINT_DEFAULT "Linting is turned off by default, use the \"lint\" target or set \ @@ -200,6 +201,26 @@ if (NOT PARTHENON_DISABLE_HDF5) install(TARGETS HDF5_C EXPORT parthenonTargets) endif() +if (PARTHENON_ENABLE_OPENPMD) +#TODO(pgrete) add logic for serial/parallel +#TODO(pgrete) add logic for internal/external build + include(FetchContent) + set(CMAKE_POLICY_DEFAULT_CMP0077 NEW) + set(openPMD_BUILD_CLI_TOOLS OFF) + set(openPMD_BUILD_EXAMPLES OFF) + set(openPMD_BUILD_TESTING OFF) + set(openPMD_BUILD_SHARED_LIBS OFF) # precedence over BUILD_SHARED_LIBS if needed + set(openPMD_INSTALL OFF) # or instead use: + # set(openPMD_INSTALL ${BUILD_SHARED_LIBS}) # only install if used as a shared library + set(openPMD_USE_PYTHON OFF) + FetchContent_Declare(openPMD + GIT_REPOSITORY "https://github.com/openPMD/openPMD-api.git" + # we need newer than the latest 0.15.2 release to support writing attriutes from a subset of ranks + GIT_TAG "1c7d7ff") # develop as of 2024-09-02 + FetchContent_MakeAvailable(openPMD) + install(TARGETS openPMD EXPORT parthenonTargets) +endif() + # Kokkos recommendatation resulting in not using default GNU extensions set(CMAKE_CXX_EXTENSIONS OFF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 37883ddcb7af..d9f25da40988 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -191,6 +191,7 @@ add_library(parthenon outputs/history.cpp outputs/io_wrapper.cpp outputs/io_wrapper.hpp + outputs/output_attr.hpp outputs/output_utils.cpp outputs/output_utils.hpp outputs/outputs.cpp @@ -203,10 +204,14 @@ add_library(parthenon outputs/parthenon_hdf5_types.hpp outputs/parthenon_xdmf.cpp outputs/parthenon_hdf5.hpp + outputs/parthenon_opmd.cpp + outputs/parthenon_opmd.hpp outputs/parthenon_xdmf.hpp outputs/restart.hpp outputs/restart_hdf5.cpp outputs/restart_hdf5.hpp + outputs/restart_opmd.cpp + outputs/restart_opmd.hpp outputs/vtk.cpp parthenon/driver.hpp @@ -326,6 +331,10 @@ if (ENABLE_HDF5) target_link_libraries(parthenon PUBLIC HDF5_C) endif() +if (PARTHENON_ENABLE_OPENPMD) + target_link_libraries(parthenon PUBLIC openPMD::openPMD) +endif() + # For Cuda with NVCC (<11.2) and C++17 Kokkos currently does not work/compile with # relaxed-constexpr, see https://github.com/kokkos/kokkos/issues/3496 # However, Parthenon heavily relies on it and there is no harm in compiling Kokkos diff --git a/src/config.hpp.in b/src/config.hpp.in index 299fe0936978..8f67ed63bc5e 100644 --- a/src/config.hpp.in +++ b/src/config.hpp.in @@ -45,6 +45,8 @@ // defne ENABLE_HDF5 or not at all #cmakedefine ENABLE_HDF5 +#cmakedefine PARTHENON_ENABLE_OPENPMD + // define PARTHENON_DISABLE_HDF5_COMPRESSION or not at all #cmakedefine PARTHENON_DISABLE_HDF5_COMPRESSION diff --git a/src/interface/params.hpp b/src/interface/params.hpp index 4b2f4a09d0cd..e99f8561d0ae 100644 --- a/src/interface/params.hpp +++ b/src/interface/params.hpp @@ -118,6 +118,12 @@ class Params { return it->second; } + const Mutability &GetMutability(const std::string &key) const { + auto const it = myMutable_.find(key); + PARTHENON_REQUIRE_THROWS(it != myMutable_.end(), "Key " + key + " doesn't exist"); + return it->second; + } + std::vector GetKeys() const { std::vector keys; for (auto &x : myParams_) { diff --git a/src/outputs/output_attr.hpp b/src/outputs/output_attr.hpp new file mode 100644 index 000000000000..8607157e09a2 --- /dev/null +++ b/src/outputs/output_attr.hpp @@ -0,0 +1,47 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2023-2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#ifndef OUTPUTS_OUTPUT_ATTR_HPP_ +#define OUTPUTS_OUTPUT_ATTR_HPP_ + +#include + +// JMM: This could probably be done with template magic but I think +// using a macro is honestly the simplest and cleanest solution here. +// Template solution would be to define a variatic class to conain the +// list of types and then a hierarchy of structs/functions to turn +// that into function calls. Preprocessor seems easier, given we're +// not manipulating this list in any way. +// The following types are the ones we allow to be stored as attributes in outputs +// (specifically within Params). +#define PARTHENON_ATTR_VALID_VEC_TYPES(T) \ + T, std::vector, ParArray1D, ParArray2D, ParArray3D, HostArray1D, \ + HostArray2D, HostArray3D, Kokkos::View, Kokkos::View, \ + ParArrayND, ParArrayHost +// JMM: This is the list of template specializations we +// "pre-instantiate" We only pre-instantiate device memory, not host +// memory. The reason is that when building with the Kokkos serial +// backend, DevMemSpace and HostMemSpace are the same and so this +// resolves to the same type in the macro, which causes problems. +#define PARTHENON_ATTR_FOREACH_VECTOR_TYPE(T) \ + PARTHENON_ATTR_APPLY(T); \ + PARTHENON_ATTR_APPLY(Kokkos::View); \ + PARTHENON_ATTR_APPLY(Kokkos::View); \ + PARTHENON_ATTR_APPLY(Kokkos::View); \ + PARTHENON_ATTR_APPLY(device_view_t) + +#endif // OUTPUTS_OUTPUT_ATTR_HPP_ diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp index e124289ac56a..3b7de06d89a0 100644 --- a/src/outputs/output_utils.cpp +++ b/src/outputs/output_utils.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -33,6 +34,7 @@ #include "mesh/meshblock.hpp" #include "outputs/output_utils.hpp" #include "parameter_input.hpp" +#include "utils/mpi_types.hpp" namespace parthenon { namespace OutputUtils { @@ -253,6 +255,47 @@ std::vector ComputeDerefinementCount(Mesh *pm) { }); } +template +std::vector FlattendedLocalToGlobal(Mesh *pm, const std::vector &data_local) { + const int n_blocks_global = pm->nbtotal; + const int n_blocks_local = static_cast(pm->block_list.size()); + + const int n_elem = data_local.size() / n_blocks_local; + PARTHENON_REQUIRE_THROWS(data_local.size() % n_blocks_local == 0, + "Results from flattened input vector does not evenly divide " + "into number of local blocks."); + std::vector data_global(n_elem * n_blocks_global); + + std::vector counts(Globals::nranks); + std::vector offsets(Globals::nranks); + + const auto &nblist = pm->GetNbList(); + counts[0] = n_elem * nblist[0]; + offsets[0] = 0; + for (int r = 1; r < Globals::nranks; r++) { + counts[r] = n_elem * nblist[r]; + offsets[r] = offsets[r - 1] + counts[r - 1]; + } + +#ifdef MPI_PARALLEL + PARTHENON_MPI_CHECK(MPI_Allgatherv(data_local.data(), counts[Globals::my_rank], + MPITypeMap::type(), data_global.data(), + counts.data(), offsets.data(), MPITypeMap::type(), + MPI_COMM_WORLD)); +#else + return data_local; +#endif + return data_global; +} + +// explicit template instantiation +template std::vector +FlattendedLocalToGlobal(Mesh *pm, const std::vector &data_local); +template std::vector +FlattendedLocalToGlobal(Mesh *pm, const std::vector &data_local); +template std::vector FlattendedLocalToGlobal(Mesh *pm, + const std::vector &data_local); + // TODO(JMM): I could make this use the other loop // functionality/high-order functions. but it was more code than this // for, I think, little benefit. @@ -329,6 +372,35 @@ std::size_t MPISum(std::size_t val) { return val; } +VariableVector GetVarsToWrite(const std::shared_ptr pmb, + const bool restart, + const std::vector &variables) { + const auto &var_vec = pmb->meshblock_data.Get()->GetVariableVector(); + auto vars_to_write = GetAnyVariables(var_vec, variables); + if (restart) { + // get all vars with flag Independent OR restart + auto restart_vars = GetAnyVariables( + var_vec, {parthenon::Metadata::Independent, parthenon::Metadata::Restart}); + for (auto restart_var : restart_vars) { + vars_to_write.emplace_back(restart_var); + } + } + return vars_to_write; +} + +std::vector GetAllVarsInfo(const VariableVector &vars, + const IndexShape &cellbounds) { + std::vector all_vars_info; + for (auto &v : vars) { + all_vars_info.emplace_back(v, cellbounds); + } + + // sort alphabetically + std::sort(all_vars_info.begin(), all_vars_info.end(), + [](const VarInfo &a, const VarInfo &b) { return a.label < b.label; }); + return all_vars_info; +} + void CheckParameterInputConsistent(ParameterInput *pin) { #ifdef MPI_PARALLEL CheckMPISizeT(); diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index e9fcee04c986..9cf01e2a0388 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -213,8 +213,8 @@ struct SwarmInfo { std::size_t count_on_rank = 0; // per-meshblock std::size_t global_offset; // global std::size_t global_count; // global - std::vector counts; // per-meshblock - std::vector offsets; // global + std::vector counts; // on local meshblocks + std::vector offsets; // global offset for local meshblocks // std::vector> masks; // used for reading swarms without defrag std::vector max_indices; // JMM: If we defrag, unneeded? void AddOffsets(const SP_Swarm &swarm); // sets above metadata @@ -236,7 +236,7 @@ struct SwarmInfo { // Copies swarmvar to host in prep for output template std::vector FillHostBuffer(const std::string vname, - ParticleVariableVector &swmvarvec) { + const ParticleVariableVector &swmvarvec) const { const auto &vinfo = var_info.at(vname); std::vector host_data(count_on_rank * vinfo.nvar); std::size_t ivec = 0; @@ -245,6 +245,7 @@ struct SwarmInfo { for (int n4 = 0; n4 < vinfo.GetN(4); ++n4) { for (int n3 = 0; n3 < vinfo.GetN(3); ++n3) { for (int n2 = 0; n2 < vinfo.GetN(2); ++n2) { + // TODO(pgrete) understand what's doing on with the blocks here... std::size_t block_idx = 0; for (auto &swmvar : swmvarvec) { // Copied extra times. JMM: If we defrag, unneeded? @@ -344,6 +345,11 @@ std::vector ComputeLocs(Mesh *pm); std::vector ComputeIDsAndFlags(Mesh *pm); std::vector ComputeDerefinementCount(Mesh *pm); +// Takes a vector containing flattened data of all rank local blocks and returns the +// flattened data over all blocks. +template +std::vector FlattendedLocalToGlobal(Mesh *pm, const std::vector &data_local); + // TODO(JMM): Potentially unsafe if MPI_UNSIGNED_LONG_LONG isn't a size_t // however I think it's probably safe to assume we'll be on systems // where this is the case? @@ -351,6 +357,17 @@ std::vector ComputeDerefinementCount(Mesh *pm); std::size_t MPIPrefixSum(std::size_t local, std::size_t &tot_count); std::size_t MPISum(std::size_t local); +// Return all variables to write, i.e., for restarts all indpendent variables and ones +// with explicit Restart flag, but also variables explicitly defined to output in the +// input file. +VariableVector GetVarsToWrite(const std::shared_ptr pmb, + const bool restart, + const std::vector &variables); + +// Returns a sorted vector of VarInfo associated with vars +std::vector GetAllVarsInfo(const VariableVector &vars, + const IndexShape &cellbounds); + void CheckParameterInputConsistent(ParameterInput *pin); } // namespace OutputUtils } // namespace parthenon diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index 430ccdf3026c..9deb80d6a2a4 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -214,8 +214,13 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { // set output variable and optional data format string used in formatted writes if ((op.file_type != "hst") && (op.file_type != "rst") && (op.file_type != "ascent") && (op.file_type != "histogram")) { - op.variables = pin->GetOrAddVector(pib->block_name, "variables", - std::vector()); + // differentiating here whether a block exists or not to not add an empty + // parameter to the input file (which might interfere with restarts) + if (pin->DoesParameterExist(pib->block_name, "variables")) { + op.variables = pin->GetVector(pib->block_name, "variables"); + } else { + op.variables = std::vector(); + } // JMM: If the requested var isn't present for a given swarm, // it is simply not output. op.swarms.clear(); // Not sure this is needed @@ -263,6 +268,20 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { pnew_type = new VTKOutput(op); } else if (op.file_type == "ascent") { pnew_type = new AscentOutput(op); + } else if (op.file_type == "openpmd") { +#ifdef PARTHENON_ENABLE_OPENPMD + const auto backend_config = + pin->GetOrAddString(op.block_name, "backend_config", "default"); + + pnew_type = new OpenPMDOutput(op, backend_config); +#else + msg << "### FATAL ERROR in Outputs constructor" << std::endl + << "Executable not configured for OpenPMD outputs, but OpenPMD file format " + << "is requested in output/restart block '" << op.block_name << "'. " + << "You can disable this block without deleting it by setting a dt < 0." + << std::endl; + PARTHENON_FAIL(msg); +#endif // ifdef PARTHENON_ENABLE_OPENPMD } else if (op.file_type == "histogram") { #ifdef ENABLE_HDF5 pnew_type = new HistogramOutput(op, pin); diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 9d37245eabc8..58f5646e4e67 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include "Kokkos_ScatterView.hpp" @@ -171,6 +172,22 @@ class AscentOutput : public OutputType { ParArray1D ghost_mask_; }; +//---------------------------------------------------------------------------------------- +//! \class OpenPMDOutput +// \brief derived OutputType class for OpenPMD based output + +class OpenPMDOutput : public OutputType { + public: + explicit OpenPMDOutput(const OutputParameters &oparams, std::string backend_config) + : OutputType(oparams), backend_config_(std::move(backend_config)) {} + void WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) override; + + private: + // path to file containing config passed to backend + std::string backend_config_; +}; + #ifdef ENABLE_HDF5 //---------------------------------------------------------------------------------------- //! \class PHDF5Output diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index cd5009d490d9..c4780e6c243d 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -360,22 +360,21 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm const auto &pmb = pm->block_list[b_idx]; bool is_allocated = false; int dealloc_count = 0; - // for each variable that this local meshblock actually has - const auto vars = get_vars(pmb); - for (auto &v : vars) { - // For reference, if we update the logic here, there's also - // a similar block in parthenon_manager.cpp - if (v->IsAllocated() && (var_name == v->label())) { - auto v_h = v->data.GetHostMirrorAndCopy(); - OutputUtils::PackOrUnpackVar( - vinfo, output_params.include_ghost_zones, index, - [&](auto index, int topo, int t, int u, int v, int k, int j, int i) { - tmpData[index] = static_cast(v_h(topo, t, u, v, k, j, i)); - }); - is_allocated = true; - dealloc_count = v->dealloc_count; - break; - } + // TODO(reviewers) Why was the loop originally there? Does the direct Get causes + // issue? + auto v = pmb->meshblock_data.Get()->GetVarPtr(var_name); + // For reference, if we update the logic here, there's also + // a similar block in parthenon_manager.cpp + if (v->IsAllocated() && (var_name == v->label())) { + auto v_h = v->data.GetHostMirrorAndCopy(); + OutputUtils::PackOrUnpackVar( + vinfo, output_params.include_ghost_zones, index, + [&](auto index, int topo, int t, int u, int v, int k, int j, int i) { + tmpData[index] = static_cast(v_h(topo, t, u, v, k, j, i)); + }); + + is_allocated = true; + dealloc_count = v->dealloc_count; } if (vinfo.is_sparse) { diff --git a/src/outputs/parthenon_hdf5.hpp b/src/outputs/parthenon_hdf5.hpp index e99c58b14815..2c0d24052da8 100644 --- a/src/outputs/parthenon_hdf5.hpp +++ b/src/outputs/parthenon_hdf5.hpp @@ -21,29 +21,9 @@ #include "defs.hpp" #include "kokkos_abstraction.hpp" +#include "output_attr.hpp" #include "parthenon_arrays.hpp" -// JMM: This could probably be done with template magic but I think -// using a macro is honestly the simplest and cleanest solution here. -// Template solution would be to define a variatic class to conain the -// list of types and then a hierarchy of structs/functions to turn -// that into function calls. Preprocessor seems easier, given we're -// not manipulating this list in any way. -#define PARTHENON_ATTR_VALID_VEC_TYPES(T) \ - T, std::vector, ParArray1D, ParArray2D, ParArray3D, HostArray1D, \ - HostArray2D, HostArray3D, Kokkos::View, Kokkos::View, \ - ParArrayND, ParArrayHost -// JMM: This is the list of template specializations we -// "pre-instantiate" We only pre-instantiate device memory, not host -// memory. The reason is that when building with the Kokkos serial -// backend, DevMemSpace and HostMemSpace are the same and so this -// resolves to the same type in the macro, which causes problems. -#define PARTHENON_ATTR_FOREACH_VECTOR_TYPE(T) \ - PARTHENON_ATTR_APPLY(T); \ - PARTHENON_ATTR_APPLY(Kokkos::View); \ - PARTHENON_ATTR_APPLY(Kokkos::View); \ - PARTHENON_ATTR_APPLY(Kokkos::View); \ - PARTHENON_ATTR_APPLY(device_view_t) // Only proceed if HDF5 output enabled #ifdef ENABLE_HDF5 diff --git a/src/outputs/parthenon_opmd.cpp b/src/outputs/parthenon_opmd.cpp new file mode 100644 index 000000000000..c0ee2aa434b2 --- /dev/null +++ b/src/outputs/parthenon_opmd.cpp @@ -0,0 +1,672 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 2024. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== +//! \file parthenon_openpmd.cpp +// \brief Output for OpenPMD https://www.openpmd.org/ (supporting various backends) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Parthenon headers +#include "basic_types.hpp" +#include "coordinates/coordinates.hpp" +#include "defs.hpp" +#include "driver/driver.hpp" +#include "globals.hpp" +#include "interface/state_descriptor.hpp" +#include "interface/swarm_default_names.hpp" +#include "interface/variable_state.hpp" +#include "mesh/mesh.hpp" +#include "mesh/meshblock.hpp" +#include "openPMD/Dataset.hpp" +#include "openPMD/Datatype.hpp" +#include "openPMD/IO/Access.hpp" +#include "openPMD/Iteration.hpp" +#include "openPMD/Mesh.hpp" +#include "openPMD/ParticleSpecies.hpp" +#include "openPMD/Series.hpp" +#include "outputs/output_attr.hpp" +#include "outputs/output_utils.hpp" +#include "outputs/outputs.hpp" +#include "outputs/parthenon_opmd.hpp" +#include "parthenon_array_generic.hpp" +#include "utils/error_checking.hpp" +#include "utils/instrument.hpp" + +// OpenPMD headers +#ifdef PARTHENON_ENABLE_OPENPMD +#include +#endif // ifdef PARTHENON_ENABLE_OPENPMD + +namespace parthenon { + +using namespace OutputUtils; + +namespace OpenPMDUtils { + +template +auto GetFlatHostVecFromView(T view) { + // Take a view and return a vector containing rank and dims and a flattened (1D) + // std::vector that can then easily be passed to OpenPMD. + // Note, this function is not + // optimial as multiple (unnecessary) copies may be done. PG didn't come up with a + // smarter way but thinks that it's not a performance issue as this is only called for + // outputs (thus not that often) and for mostly small amounts of data. With a C++20 span + // we could probably direct reuse the host mirror data pointer. + auto view_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), view); + + using base_t = typename std::remove_pointer::type; + auto host_vec = std::vector(view_h.size()); + for (auto i = 0; i < view_h.size(); i++) { + host_vec[i] = view_h.data()[i]; + } + // cpplint demands compile constants be all caps + constexpr auto RANK = static_cast(T::rank); + std::vector rank_and_dims(RANK + 1); + rank_and_dims[0] = RANK; + for (size_t d = 0; d < RANK; ++d) { + rank_and_dims[1 + d] = view.extent_int(d); + } + return std::make_tuple(rank_and_dims, host_vec); +} + +template +void WriteAllParamsOfType(const Params ¶ms, const std::string &prefix, + openPMD::Iteration *it) { + for (const auto &key : params.GetKeys()) { + const auto type = params.GetType(key); + if (type == std::type_index(typeid(T))) { + auto full_path = prefix + delim + key; + // The '/' is kind of a reserved character in the OpenPMD standard, which results + // in attribute keys with said character not being exposed. + // Thus we replace it. + std::replace(full_path.begin(), full_path.end(), '/', delim[0]); + + if constexpr (implements::value) { + const auto &view = params.Get(key); + auto [rank_and_dims, host_vec] = GetFlatHostVecFromView(view); + it->setAttribute(full_path + ".rankdims", rank_and_dims); + it->setAttribute(full_path, host_vec); + } else if constexpr (is_specialization_of::value) { + const auto &view = params.Get(key).KokkosView(); + auto [rank_and_dims, host_vec] = GetFlatHostVecFromView(view); + it->setAttribute(full_path + ".rankdims", rank_and_dims); + it->setAttribute(full_path, host_vec); + } else { + it->setAttribute(full_path, params.Get(key)); + } + } + } +} + +template +void WriteAllParamsOfMultipleTypes(const Params ¶ms, const std::string &prefix, + openPMD::Iteration *it) { + ([&] { WriteAllParamsOfType(params, prefix, it); }(), ...); +} + +template +void WriteAllParams(const Params ¶ms, const std::string &prefix, + openPMD::Iteration *it) { + WriteAllParamsOfMultipleTypes(params, prefix, it); +} + +void WriteAllParams(const Params ¶ms, const std::string &pkg_name, + openPMD::Iteration *it) { + using OpenPMDUtils::delim; + const std::string prefix = "Params" + delim + pkg_name; + // check why this (vector of bool) doesn't work + // WriteAllParams(params, prefix, it); + WriteAllParamsOfType(params, prefix, it); + WriteAllParams(params, prefix, it); + WriteAllParams(params, prefix, it); + WriteAllParams(params, prefix, it); + WriteAllParams(params, prefix, it); + WriteAllParams(params, prefix, it); + WriteAllParams(params, prefix, it); + WriteAllParams(params, prefix, it); +} + +template +void WriteSwarmVar(const SwarmInfo &swinfo, openPMD::ParticleSpecies swm, + openPMD::Iteration it) { + auto &vars_of_type_T = std::get>(swinfo.vars); + for (const auto &[vname, swmvarvec] : vars_of_type_T) { + const auto &vinfo = swinfo.var_info.at(vname); + auto host_data = swinfo.FillHostBuffer(vname, swmvarvec); + + auto const dataset = openPMD::Dataset(openPMD::determineDatatype(host_data.data()), + {swinfo.global_count}); + // TODO(pgrete) ask OpenPMD group if this is the right approach of if our non-scalar + // partices should be a multi-D `dataset` if is scalar + if (vinfo.tensor_rank == 0) { + // special sauce to align "positions" with standard + std::string particle_record; + std::string particle_record_component; + if (vname == swarm_position::x::name()) { + particle_record = "position"; + particle_record_component = "x"; + } else if (vname == swarm_position::y::name()) { + particle_record = "position"; + particle_record_component = "y"; + } else if (vname == swarm_position::z::name()) { + particle_record = "position"; + particle_record_component = "z"; + } else { + particle_record = vname; + particle_record_component = openPMD::MeshRecordComponent::SCALAR; + } + openPMD::RecordComponent rc = swm[particle_record][particle_record_component]; + rc.resetDataset(dataset); + // only write if there's sth to write (otherwise the host_data nullptr is caught) + if (swinfo.count_on_rank != 0) { + rc.storeChunk(host_data, {swinfo.global_offset}, {host_data.size()}); + } + + // if positional, add offsets + if (particle_record_component != openPMD::MeshRecordComponent::SCALAR) { + auto rc_offset = swm["positionOffset"][particle_record_component]; + rc_offset.resetDataset(dataset); + rc_offset.makeConstant(0.0); + } + + // else flatten components + } else { + for (auto n = 0; n < vinfo.nvar; n++) { + openPMD::RecordComponent rc = swm[vname][std::to_string(n)]; + rc.resetDataset(dataset); + // only write if there's sth to write (otherwise the host_data nullptr is caught) + if (swinfo.count_on_rank != 0) { + rc.storeChunkRaw(&host_data[n * swinfo.count_on_rank], {swinfo.global_offset}, + {swinfo.count_on_rank}); + } + } + } + // Flush because the host buffer is temporary + it.seriesFlush(); + } +} +std::tuple GetMeshRecordAndComponentNames(const VarInfo &vinfo, + const int comp_idx, + const int level) { + std::string comp_name; + if (vinfo.is_vector) { + if (comp_idx == 0) { + comp_name = "x"; + } else if (comp_idx == 1) { + comp_name = "y"; + } else if (comp_idx == 2) { + comp_name = "z"; + } else { + PARTHENON_THROW("Expected component index doesn't match vector expectation."); + } + // Current unclear how to properly handle other vectors and tensors, so everything + // that not's a proper vector is a a scalar for now. + } else { + comp_name = openPMD::MeshRecordComponent::SCALAR; + } + // TODO(pgrete) need to make sure that var names are allowed within standard + const std::string &mesh_record_name = vinfo.label + "_" + + vinfo.component_labels[comp_idx] + "_lvl" + + std::to_string(level); + // return std::make_tuple(mesh_record_name, comp_name); + return {mesh_record_name, comp_name}; +} + +std::tuple +GetChunkOffsetAndExtent(Mesh *pm, std::shared_ptr pmb) { + openPMD::Offset chunk_offset; + openPMD::Extent chunk_extent; + const auto loc = pm->Forest().GetLegacyTreeLocation(pmb->loc); + if (pm->ndim == 3) { + chunk_offset = {loc.lx3() * static_cast(pmb->block_size.nx(X3DIR)), + loc.lx2() * static_cast(pmb->block_size.nx(X2DIR)), + loc.lx1() * static_cast(pmb->block_size.nx(X1DIR))}; + chunk_extent = {static_cast(pmb->block_size.nx(X3DIR)), + static_cast(pmb->block_size.nx(X2DIR)), + static_cast(pmb->block_size.nx(X1DIR))}; + } else if (pm->ndim == 2) { + chunk_offset = {loc.lx2() * static_cast(pmb->block_size.nx(X2DIR)), + loc.lx1() * static_cast(pmb->block_size.nx(X1DIR))}; + chunk_extent = {static_cast(pmb->block_size.nx(X2DIR)), + static_cast(pmb->block_size.nx(X1DIR))}; + } else { + PARTHENON_THROW("1D output for openpmd not yet supported."); + } + return {chunk_offset, chunk_extent}; +} +} // namespace OpenPMDUtils + +//---------------------------------------------------------------------------------------- +//! \fn void OpenPMDOutput:::WriteOutputFile(Mesh *pm) +// \brief Expose mesh and all Cell variables for processing with Ascent +void OpenPMDOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) { +#ifndef PARTHENON_ENABLE_OPENPMD + if (Globals::my_rank == 0) { + PARTHENON_WARN("OpenPMD output requested by input file, but OpenPMD support not " + "compiled in. Skipping this output type."); + } +#else + using openPMD::Access; + using openPMD::Series; + + // TODO(pgrete) .h5 for hd5 and .bp for ADIOS2 or .json for JSON + // TODO(pgrete) check if CREATE is the correct pattern (for not overwriting the series + // but an interation) This just describes the pattern of the filename. The correct file + // will be accessed through the iteration idx below. The file suffix maps to the chosen + // backend. + // TODO(pgrete) add final and now logic + // Prepending @ indicates that the config is a file to be read and parsed. + std::string backend_config = + backend_config_ == "default" ? "{}" : "@" + backend_config_; + + auto filename = output_params.file_basename + "." + output_params.file_id; + if (signal == SignalHandler::OutputSignal::now) { + filename.append(".now"); + } + filename.append(".%05T.bp"); + Series series = Series(filename, Access::CREATE, +#ifdef MPI_PARALLEL + MPI_COMM_WORLD, +#endif + backend_config); + + // TODO(pgrete) How to handle downstream info, e.g., on how/what defines a vector? + // TODO(pgrete) Should we update for restart or only set this once? Or make it per + // iteration? + // ... = pin->GetString(output_params.block_name, "actions_file"); + series.setAuthor("My Name file naming + auto it = series.iterations[output_params.file_number]; + it.open(); // explicit open() is important when run in parallel + + auto const &first_block = *(pm->block_list.front()); + + // TODO(?) in principle, we could abstract this to a more general WriteAttributes place + // and reuse for hdf5 and OpenPMD output with corresponing calls + // -------------------------------------------------------------------------------- // + // WRITING ATTRIBUTES // + // -------------------------------------------------------------------------------- // + + // Note, that profiling is likely skewed as data is actually written to disk/flushed + // only later. + Kokkos::Profiling::pushRegion("write Attributes"); + // First the ones required by the OpenPMD standard + if (tm != nullptr) { + it.setTime(tm->time); + it.setDt(tm->dt); + it.setAttribute("NCycle", tm->ncycle); + } else { + it.setTime(-1.0); + it.setDt(-1.0); + } + { + PARTHENON_INSTRUMENT_REGION("Dump Params"); + + for (const auto &[pkg_name, pkg] : pm->packages.AllPackages()) { + const auto ¶ms = pkg->AllParams(); + OpenPMDUtils::WriteAllParams(params, pkg_name, &it); + } + } + // Then our own + { + PARTHENON_INSTRUMENT_REGION("write input"); + // write input key-value pairs + std::ostringstream oss; + pin->ParameterDump(oss); + it.setAttribute("InputFile", oss.str()); + } + + { + // It's not clear we need all these attributes, but they mirror what's done in the + // hdf5 output. + it.setAttribute("WallTime", Driver::elapsed_main()); + it.setAttribute("NumDims", pm->ndim); + it.setAttribute("NumMeshBlocks", pm->nbtotal); + it.setAttribute("MaxLevel", pm->GetCurrentLevel() - pm->GetRootLevel()); + // write whether we include ghost cells or not + it.setAttribute("IncludesGhost", output_params.include_ghost_zones ? 1 : 0); + // write number of ghost cells in simulation + it.setAttribute("NGhost", Globals::nghost); + it.setAttribute("Coordinates", std::string(first_block.coords.Name()).c_str()); + + // restart info, write always + it.setAttribute("NBNew", pm->nbnew); + it.setAttribute("NBDel", pm->nbdel); + it.setAttribute("RootLevel", pm->GetLegacyTreeRootLevel()); + it.setAttribute("Refine", pm->adaptive ? 1 : 0); + it.setAttribute("Multilevel", pm->multilevel ? 1 : 0); + + it.setAttribute("BlocksPerPE", pm->GetNbList()); + + // Mesh block size + const auto base_block_size = pm->GetDefaultBlockSize(); + it.setAttribute("MeshBlockSize", + std::vector{base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), + base_block_size.nx(X3DIR)}); + + // RootGridDomain - float[9] array with xyz mins, maxs, rats (dx(i)/dx(i-1)) + it.setAttribute( + "RootGridDomain", + std::vector{pm->mesh_size.xmin(X1DIR), pm->mesh_size.xmax(X1DIR), + pm->mesh_size.xrat(X1DIR), pm->mesh_size.xmin(X2DIR), + pm->mesh_size.xmax(X2DIR), pm->mesh_size.xrat(X2DIR), + pm->mesh_size.xmin(X3DIR), pm->mesh_size.xmax(X3DIR), + pm->mesh_size.xrat(X3DIR)}); + + // Root grid size (number of cells at root level) + it.setAttribute("RootGridSize", + std::vector{pm->mesh_size.nx(X1DIR), pm->mesh_size.nx(X2DIR), + pm->mesh_size.nx(X3DIR)}); + + // Boundary conditions + auto arr_to_vec = [](const auto &arr) { + std::vector vec(BOUNDARY_NFACES); + for (int i = 0; i < BOUNDARY_NFACES; i++) { + vec[i] = arr.at(i); + } + return vec; + }; + it.setAttribute("BoundaryConditions", arr_to_vec(pm->mesh_bc_names)); + it.setAttribute("SwarmBoundaryConditions", arr_to_vec(pm->mesh_swarm_bc_names)); + } // Info section + + Kokkos::Profiling::popRegion(); // write Attributes + + // Write block metadata + { + // Manually gather all block data first as it allows to use the (simpler) + // Attribute interface rather than writing a distributed dataset -- especially as all + // data is being read on restart by every rank anyway. + std::vector loc_local = OutputUtils::ComputeLocs(pm); + auto loc_global = FlattendedLocalToGlobal(pm, loc_local); + it.setAttribute("loc.lx123", loc_global); + + std::vector id_local = OutputUtils::ComputeIDsAndFlags(pm); + auto id_global = FlattendedLocalToGlobal(pm, id_local); + it.setAttribute("loc.level-gid-lid-cnghost-gflag", id_global); + + // derefinement count + std::vector derefcnt_local = OutputUtils::ComputeDerefinementCount(pm); + auto derefcnt_global = FlattendedLocalToGlobal(pm, derefcnt_local); + it.setAttribute("derefinement_count", derefcnt_global); + } + + // TODO(pgrete) check var name standard compatiblity + // e.g., description: names of records and their components are only allowed to contain + // the characters a-Z, the numbers 0-9 and the underscore _ + + const int num_blocks_local = static_cast(pm->block_list.size()); + + // -------------------------------------------------------------------------------- // + // WRITING VARIABLES DATA // + // -------------------------------------------------------------------------------- // + Kokkos::Profiling::pushRegion("write all variable data"); + + auto &bounds = pm->block_list.front()->cellbounds; + // get list of all vars, just use the first block since the list is the same for all + // blocks + // TODO(pgrete) add restart_ var to output + // TODO(pgrete) check if this needs to be updated/unifed with get_var logic in hdf5 + auto all_vars_info = GetAllVarsInfo( + GetVarsToWrite(pm->block_list.front(), true, output_params.variables), bounds); + + // We're currently writing (flushing) one var at a time. This saves host memory but + // results more smaller write. Might be updated in the future. + // Allocate space for largest size variable + int var_size_max = 0; + for (auto &vinfo : all_vars_info) { + const auto var_size = vinfo.Size(); + var_size_max = std::max(var_size_max, var_size); + } + + // TODO(pgrete) adjust for single prec output + // openPMD::Datatype dtype = openPMD::determineDatatype(); + using OutT = + Real; // typename std::conditional::type; + std::vector tmp_data(var_size_max * num_blocks_local); + + // TODO(pgrete) This needs to be in the loop for non-cell-centered vars + auto ib = bounds.GetBoundsI(IndexDomain::interior); + auto jb = bounds.GetBoundsJ(IndexDomain::interior); + auto kb = bounds.GetBoundsK(IndexDomain::interior); + // for each variable we write + for (auto &vinfo : all_vars_info) { + PARTHENON_INSTRUMENT_REGION("Write variable loop") + + // Reset host write bufer. Not really necessary, but doesn't hurt. + memset(tmp_data.data(), 0, tmp_data.size() * sizeof(OutT)); + uint64_t tmp_offset = 0; + + if (vinfo.is_vector) { + // sanity check + PARTHENON_REQUIRE_THROWS( + vinfo.GetDim(4) == pm->ndim && vinfo.GetDim(5) == 1 && vinfo.GetDim(6) == 1, + "A 'standard' vector is expected to only have components matching the " + "dimensionality of the simulation.") + } + + for (auto &pmb : pm->block_list) { + // TODO(pgrete) check if we should skip the suffix for level 0 + const auto level = pmb->loc.level() - pm->GetRootLevel(); + + for (int comp_idx = 0; comp_idx < vinfo.component_labels.size(); comp_idx++) { + const auto [record_name, comp_name] = + OpenPMDUtils::GetMeshRecordAndComponentNames(vinfo, comp_idx, level); + + // Create the mesh_record for this variable at the given level (if it doesn't + // exist yet) + if (!it.meshes.contains(record_name)) { + auto mesh_record = it.meshes[record_name]; + + // These following attributes are shared across all components of the record. + + PARTHENON_REQUIRE_THROWS( + typeid(Coordinates_t) == typeid(UniformCartesian), + "OpenPMD in Parthenon currently only supports Cartesian coordinates."); + mesh_record.setGeometry(openPMD::Mesh::Geometry::cartesian); + auto &coords = pmb->coords; + // For uniform Cartesian, all dxN are const across the block so we just pick the + // first index. + Real dx1 = coords.CellWidth(0, 0, 0); + Real dx2 = coords.CellWidth(0, 0, 0); + Real dx3 = coords.CellWidth(0, 0, 0); + + // TODO(pgrete) check if this should be tied to the MemoryLayout + mesh_record.setDataOrder(openPMD::Mesh::DataOrder::C); + + auto mesh_comp = mesh_record[comp_name]; + // TODO(pgrete) needs to be updated for face and edges etc + // Also this feels wrong for deep hierachies... + auto effective_nx = static_cast(std::pow(2, level)); + openPMD::Extent global_extent; + if (pm->ndim == 3) { + mesh_record.setGridSpacing(std::vector{dx3, dx2, dx1}); + mesh_record.setAxisLabels({"z", "y", "x"}); + mesh_record.setGridGlobalOffset({ + pm->mesh_size.xmin(X3DIR), + pm->mesh_size.xmin(X2DIR), + pm->mesh_size.xmin(X1DIR), + }); + // TODO(pgrete) needs to be updated for face and edges etc + mesh_comp.setPosition(std::vector{0.5, 0.5, 0.5}); + global_extent = { + static_cast(pm->mesh_size.nx(X3DIR)) * effective_nx, + static_cast(pm->mesh_size.nx(X2DIR)) * effective_nx, + static_cast(pm->mesh_size.nx(X1DIR)) * effective_nx, + }; + } else if (pm->ndim == 2) { + mesh_record.setGridSpacing(std::vector{dx2, dx1}); + mesh_record.setAxisLabels({"y", "x"}); + mesh_record.setGridGlobalOffset({ + pm->mesh_size.xmin(X2DIR), + pm->mesh_size.xmin(X1DIR), + }); + + // TODO(pgrete) needs to be updated for face and edges etc + mesh_comp.setPosition(std::vector{0.5, 0.5}); + global_extent = { + static_cast(pm->mesh_size.nx(X2DIR)) * effective_nx, + static_cast(pm->mesh_size.nx(X1DIR)) * effective_nx, + }; + + } else { + PARTHENON_THROW("1D output for openpmd not yet supported."); + } + // Handling this here to now re-reset dataset later when iterating through the + // blocks + auto const dataset = + openPMD::Dataset(openPMD::determineDatatype(), global_extent); + // TODO(pgrete) check whether this should/need to be a collective so that the + // mesh generation should be done across all ranks prior to writing data, rather + // than in-situ for the local blocks only + mesh_comp.resetDataset(dataset); + + // TODO(pgrete) need unitDimension and timeOffset for this record? + } + } + + // Now that the mesh record exists, actually write the data + auto out_var = pmb->meshblock_data.Get()->GetVarPtr(vinfo.label); + PARTHENON_REQUIRE_THROWS(out_var->metadata().Where() == + MetadataFlag(Metadata::Cell), + "Currently only cell centered vars are supported."); + + if (out_var->IsAllocated()) { + // TODO(pgrete) check if we can work with a direct copy from a subview to not + // duplicate the memory footprint here +#if 0 + // Pick a subview of the active cells of this component + auto const data = Kokkos::subview( + var->data, 0, 0, icomp, std::make_pair(kb.s, kb.e + 1), + std::make_pair(jb.s, jb.e + 1), std::make_pair(ib.s, ib.e + 1)); + + // Map a view onto a host allocation (so that we can call deep_copy) + auto component_buffer = buffer_list.emplace_back(ncells); + Kokkos::View> + component_buffer_view(component_buffer.data(), nk, nj, ni); + Kokkos::deep_copy(component_buffer_view, data); +#endif + auto out_var_h = out_var->data.GetHostMirrorAndCopy(); + int comp_idx = 0; + const auto &Nt = out_var->GetDim(6); + const auto &Nu = out_var->GetDim(5); + const auto &Nv = out_var->GetDim(4); + // loop over all components + for (int t = 0; t < Nt; ++t) { + for (int u = 0; u < Nu; ++u) { + for (int v = 0; v < Nv; ++v) { + const auto [record_name, comp_name] = + OpenPMDUtils::GetMeshRecordAndComponentNames(vinfo, comp_idx, level); + auto mesh_comp = it.meshes[record_name][comp_name]; + + const auto comp_offset = tmp_offset; + for (int k = kb.s; k <= kb.e; ++k) { + for (int j = jb.s; j <= jb.e; ++j) { + for (int i = ib.s; i <= ib.e; ++i) { + tmp_data[tmp_offset] = static_cast(out_var_h(t, u, v, k, j, i)); + tmp_offset++; + } + } + } + const auto [chunk_offset, chunk_extent] = + OpenPMDUtils::GetChunkOffsetAndExtent(pm, pmb); + mesh_comp.storeChunkRaw(&tmp_data[comp_offset], chunk_offset, chunk_extent); + comp_idx += 1; + } + } + } // loop over components + } // out_var->IsAllocated() + } // loop over blocks + it.seriesFlush(); + } // loop over vars + Kokkos::Profiling::popRegion(); // write all variable data + + // -------------------------------------------------------------------------------- // + // WRITING PARTICLE DATA // + // -------------------------------------------------------------------------------- // + + Kokkos::Profiling::pushRegion("write particle data"); + // TODO(pgrete) as above, first wrt differentiating between restart_ (last arg) + AllSwarmInfo all_swarm_info(pm->block_list, output_params.swarms, true); + for (auto &[swname, swinfo] : all_swarm_info.all_info) { + openPMD::ParticleSpecies swm = it.particles[swname]; + // These indicate particles/meshblock and location in global index + // space where each meshblock starts + auto counts_global = FlattendedLocalToGlobal(pm, swinfo.counts); + swm.setAttribute("counts", counts_global); + auto offsets_global = FlattendedLocalToGlobal(pm, swinfo.offsets); + swm.setAttribute("offsets", offsets_global); + + if (swinfo.global_count == 0) { + continue; + } + + OpenPMDUtils::WriteSwarmVar(swinfo, swm, it); + OpenPMDUtils::WriteSwarmVar(swinfo, swm, it); + + // From the HDF5 output: + // If swarm does not contain an "id" object, generate a sequential + // one for vis. + // BUT PG: this may break things in unpredicable ways + // I'm in favor of enforcing a global id somehow. We shold discuss. + PARTHENON_REQUIRE_THROWS(swinfo.var_info.count("id") != 0, + "Particles should always carry a unique, persistent id!"); + } + Kokkos::Profiling::popRegion(); // write particle data + + // The iteration can be closed in order to help free up resources. + // The iteration's content will be flushed automatically. + // An iteration once closed cannot (yet) be reopened. + it.close(); + series.close(); +#endif // ifndef PARTHENON_ENABLE_OPENPMD + + // advance output parameters if this is not a triggered (now or final) output + if (signal == SignalHandler::OutputSignal::none) { + output_params.file_number++; + output_params.next_time += output_params.dt; + pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); + pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + } +} + +} // namespace parthenon diff --git a/src/outputs/parthenon_opmd.hpp b/src/outputs/parthenon_opmd.hpp new file mode 100644 index 000000000000..9e0c04e99a4c --- /dev/null +++ b/src/outputs/parthenon_opmd.hpp @@ -0,0 +1,51 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +#ifndef OUTPUTS_PARTHENON_OPMD_HPP_ +#define OUTPUTS_PARTHENON_OPMD_HPP_ +//! \file restart_opmd.hpp +// \brief Provides support for restarting from OpenPMD output + +#include +#include +#include + +#include "mesh/meshblock.hpp" +#include "openPMD/Dataset.hpp" +#include "openPMD/Iteration.hpp" +#include "outputs/output_utils.hpp" + +namespace parthenon { + +namespace OpenPMDUtils { + +template +void RestoreViewAttribute(const std::string &full_path, T &view, openPMD::Iteration *it); + +void WriteAllParams(const Params ¶ms, const std::string &prefix, + openPMD::Iteration *it); + +// Deliminter to separate packages and parameters in attributes. +// More or less a workaround as the OpenPMD API does currently not expose +// access to non-standard groups (such as "Params" versus the standard "meshes"). +inline static const std::string delim = "~"; + +// Construct OpenPMD Mesh "record" name and comonnent identifier. +// - comp_idx is a flattended index over all components of the vectors and tensors, i.e., +// the typical v,u,t indices. +// - level is the current effective level of the Mesh record +std::tuple +GetMeshRecordAndComponentNames(const OutputUtils::VarInfo &vinfo, const int comp_idx, + const int level); + +// Calculate logical location on effective mesh (i.e., a mesh with size that matches full +// coverage at given resolution on a particular level) +// TODO(pgrete) needs to be updated to properly work with Forests +std::tuple +GetChunkOffsetAndExtent(Mesh *pm, std::shared_ptr pmb); + +} // namespace OpenPMDUtils +} // namespace parthenon +#endif // OUTPUTS_PARTHENON_OPMD_HPP_ diff --git a/src/outputs/restart.hpp b/src/outputs/restart.hpp index 138949dd1667..3784f4a2dc0d 100644 --- a/src/outputs/restart.hpp +++ b/src/outputs/restart.hpp @@ -106,7 +106,7 @@ class RestartReader { // fills internal data for given pointer virtual void ReadBlocks(const std::string &name, IndexRange range, const OutputUtils::VarInfo &info, std::vector &dataVec, - int file_output_format_version) const = 0; + int file_output_format_version, Mesh *pmesh) const = 0; // Gets the data from a swarm var on current rank. Assumes all // blocks are contiguous. Fills dataVec based on shape from swarmvar diff --git a/src/outputs/restart_hdf5.cpp b/src/outputs/restart_hdf5.cpp index 8078fbd5fa9a..3c1396466ec6 100644 --- a/src/outputs/restart_hdf5.cpp +++ b/src/outputs/restart_hdf5.cpp @@ -207,7 +207,8 @@ void RestartReaderHDF5::ReadParams(const std::string &name, Params &p) { void RestartReaderHDF5::ReadBlocks(const std::string &name, IndexRange range, const OutputUtils::VarInfo &info, std::vector &dataVec, - int file_output_format_version) const { + int file_output_format_version, + Mesh * /*pmesh*/) const { #ifndef ENABLE_HDF5 PARTHENON_FAIL("Restart functionality is not available because HDF5 is disabled"); #else // HDF5 enabled @@ -246,6 +247,7 @@ void RestartReaderHDF5::ReadBlocks(const std::string &name, IndexRange range, std::to_string(total_count) + ")"); const H5S memspace = H5S::FromHIDCheck(H5Screate_simple(total_dim, count, NULL)); + // TODO(reviewer) What's going on here? The follow line is identical to the one above. PARTHENON_HDF5_CHECK( H5Sselect_hyperslab(hdl.dataspace, H5S_SELECT_SET, offset, NULL, count, NULL)); diff --git a/src/outputs/restart_hdf5.hpp b/src/outputs/restart_hdf5.hpp index 4c346fdac78d..41e11d26787e 100644 --- a/src/outputs/restart_hdf5.hpp +++ b/src/outputs/restart_hdf5.hpp @@ -115,7 +115,7 @@ class RestartReaderHDF5 : public RestartReader { // fills internal data for given pointer void ReadBlocks(const std::string &name, IndexRange range, const OutputUtils::VarInfo &info, std::vector &dataVec, - int file_output_format_version) const override; + int file_output_format_version, Mesh *pmesh) const override; // Gets the data from a swarm var on current rank. Assumes all // blocks are contiguous. Fills dataVec based on shape from swarmvar diff --git a/src/outputs/restart_opmd.cpp b/src/outputs/restart_opmd.cpp new file mode 100644 index 000000000000..11f9943124a1 --- /dev/null +++ b/src/outputs/restart_opmd.cpp @@ -0,0 +1,231 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file restart_opmd.cpp +// \brief Restarts a simulation from an OpenPMD output with ADIOS2 backend + +#include +#include +#include +#include +#include +#include +#include + +#include "basic_types.hpp" +#include "interface/params.hpp" +#include "openPMD/Iteration.hpp" +#include "openPMD/Series.hpp" +#include "outputs/output_attr.hpp" +#include "outputs/parthenon_opmd.hpp" +#include "outputs/restart.hpp" +#include "outputs/restart_opmd.hpp" +#include "utils/error_checking.hpp" + +namespace parthenon { + +//---------------------------------------------------------------------------------------- +//! \fn void RestartReader::RestartReader(const std::string filename) +// \brief Opens the restart file and stores appropriate file handle in fh_ +RestartReaderOPMD::RestartReaderOPMD(const char *filename) : filename_(filename) { + // This silly logic is required as the unit tests may or may not define MPI_PARALLEL but + // are always run in serial. +#ifdef MPI_PARALLEL + int mpi_initialized; + PARTHENON_MPI_CHECK(MPI_Initialized(&mpi_initialized)); + if (mpi_initialized) { + series = openPMD::Series(filename, openPMD::Access::READ_ONLY, MPI_COMM_WORLD); + } else { + series = openPMD::Series(filename, openPMD::Access::READ_ONLY); + } +#else + series = openPMD::Series(filename, openPMD::Access::READ_ONLY); + +#endif + PARTHENON_REQUIRE_THROWS( + series.iterations.size() == 1, + "Parthenon restarts should only contain one iteration/timestep."); + std::uint64_t idx; + for (const auto &i : series.iterations) { + idx = i.first; + } + it = std::make_unique(series.iterations[idx]); + // Explicitly open (important for parallel execution) + it->open(); +} + +int RestartReaderOPMD::GetOutputFormatVersion() const { + // TODO(pgrete) move info to shared header and introduce constexpr var + if (it->containsAttribute("OutputFormatVersion")) { + return it->getAttribute("OutputFormatVersion").get(); + } else { + return -1; + } +} + +RestartReaderOPMD::SparseInfo RestartReaderOPMD::GetSparseInfo() const { + // TODO(pgrete) needs impl + return {}; +} + +RestartReaderOPMD::MeshInfo RestartReaderOPMD::GetMeshInfo() const { + RestartReaderOPMD::MeshInfo mesh_info; + mesh_info.nbnew = it->getAttribute("NBNew").get(); + mesh_info.nbdel = it->getAttribute("NBDel").get(); + mesh_info.nbtotal = it->getAttribute("NumMeshBlocks").get(); + mesh_info.root_level = it->getAttribute("RootLevel").get(); + + mesh_info.block_size = it->getAttribute("MeshBlockSize").get>(); + mesh_info.includes_ghost = it->getAttribute("IncludesGhost").get(); + mesh_info.n_ghost = it->getAttribute("NGhost").get(); + + mesh_info.grid_dim = it->getAttribute("RootGridDomain").get>(); + mesh_info.lx123 = it->getAttribute("loc.lx123").get>(); + mesh_info.level_gid_lid_cnghost_gflag = + it->getAttribute("loc.level-gid-lid-cnghost-gflag").get>(); + mesh_info.derefinement_count = + it->getAttribute("derefinement_count").get>(); + + return mesh_info; +} + +SimTime RestartReaderOPMD::GetTimeInfo() const { + SimTime time_info{}; + + time_info.time = it->time(); + time_info.dt = it->dt(); + time_info.ncycle = it->getAttribute("NCycle").get(); + + return time_info; +} +// Gets the counts and offsets for MPI ranks for the meshblocks set +// by the indexrange. Returns the total count on this rank. +std::size_t RestartReaderOPMD::GetSwarmCounts(const std::string &swarm, + const IndexRange &range, + std::vector &counts, + std::vector &offsets) { + // datasets + auto counts_dset = + it->particles[swarm].getAttribute("counts").get>(); + auto offsets_dset = + it->particles[swarm].getAttribute("offsets").get>(); + + // Read data for requested blocks in range + counts.resize(range.e - range.s + 1); + offsets.resize(range.e - range.s + 1); + + std::copy(counts_dset.begin() + range.s, counts_dset.begin() + range.e + 1, + counts.begin()); + std::copy(offsets_dset.begin() + range.s, offsets_dset.begin() + range.e + 1, + offsets.begin()); + + // Compute total count rank + std::size_t total_count_on_rank = std::accumulate(counts.begin(), counts.end(), 0); + return total_count_on_rank; +} + +template +void RestartReaderOPMD::ReadAllParamsOfType(const std::string &prefix, Params ¶ms) { + for (const auto &key : params.GetKeys()) { + using OpenPMDUtils::delim; + const auto type = params.GetType(key); + auto mutability = params.GetMutability(key); + if (type == std::type_index(typeid(T)) && mutability == Params::Mutability::Restart) { + auto full_path = prefix + delim + key; + // The '/' is kind of a reserved character in the OpenPMD standard, which results + // in attribute keys with said character not being exposed. + // Thus we replace it. + std::replace(full_path.begin(), full_path.end(), '/', delim[0]); + + T val; + if constexpr (implements::value) { + val = params.Get(key); + RestoreViewAttribute(full_path, val); + } else if constexpr (is_specialization_of::value) { + val = params.Get(key); + auto &view = val.KokkosView(); + RestoreViewAttribute(full_path, view); + } else { + val = it->getAttribute(full_path).get(); + } + params.Update(key, val); + } + } +} + +template +void RestartReaderOPMD::ReadAllParamsOfMultipleTypes(const std::string &prefix, + Params &p) { + ([&] { ReadAllParamsOfType(prefix, p); }(), ...); +} + +template +void RestartReaderOPMD::ReadAllParams(const std::string &pkg_name, Params &p) { + ReadAllParamsOfMultipleTypes(pkg_name, p); +} +void RestartReaderOPMD::ReadParams(const std::string &pkg_name, Params &p) { + using OpenPMDUtils::delim; + const auto prefix = "Params" + delim + pkg_name; + ReadAllParams(prefix, p); + ReadAllParams(prefix, p); + ReadAllParams(prefix, p); + ReadAllParams(prefix, p); + ReadAllParams(prefix, p); + ReadAllParams(prefix, p); + ReadAllParams(prefix, p); + // TODO(pgrete) same as for the writing. fix vec of bool + ReadAllParamsOfType(prefix, p); +} + +void RestartReaderOPMD::ReadBlocks(const std::string &var_name, IndexRange block_range, + const OutputUtils::VarInfo &vinfo, + std::vector &data_vec, + int file_output_format_version, Mesh *pm) const { + int64_t comp_offset = 0; // offset data_vector to store component data + for (auto &pmb : pm->block_list) { + // TODO(pgrete) check if we should skip the suffix for level 0 + const auto level = pmb->loc.level() - pm->GetRootLevel(); + + int comp_idx = 0; // used in label for non-vector variables + const auto &Nt = vinfo.GetDim(6); + const auto &Nu = vinfo.GetDim(5); + const auto &Nv = vinfo.GetDim(4); + // loop over all components + for (int t = 0; t < Nt; ++t) { + for (int u = 0; u < Nu; ++u) { + for (int v = 0; v < Nv; ++v) { + // Get the correct record + const auto [record_name, comp_name] = + OpenPMDUtils::GetMeshRecordAndComponentNames(vinfo, comp_idx, level); + + PARTHENON_REQUIRE_THROWS(it->meshes.contains(record_name), + "Missing mesh record '" + record_name + + "' in restart file."); + auto mesh_record = it->meshes[record_name]; + PARTHENON_REQUIRE_THROWS(mesh_record.contains(comp_name), + "Missing component'" + comp_name + + "' in mesh record '" + record_name + + "' of restart file."); + auto mesh_comp = mesh_record[comp_name]; + + const auto [chunk_offset, chunk_extent] = + OpenPMDUtils::GetChunkOffsetAndExtent(pm, pmb); + mesh_comp.loadChunkRaw(&data_vec[comp_offset], chunk_offset, chunk_extent); + // TODO(pgrete) check if output utils machinery can be used for non-cell + // centered fields, which might not be that straightforward as a global mesh + // is stored rather than individual blocks. + comp_offset += pmb->block_size.nx(X1DIR) * pmb->block_size.nx(X2DIR) * + pmb->block_size.nx(X3DIR); + comp_idx += 1; + } + } + } // loop over components + } // loop over blocks + + // Now actually read the registered chunks form disk + it->seriesFlush(); +} + +} // namespace parthenon diff --git a/src/outputs/restart_opmd.hpp b/src/outputs/restart_opmd.hpp new file mode 100644 index 000000000000..5f79a7662012 --- /dev/null +++ b/src/outputs/restart_opmd.hpp @@ -0,0 +1,164 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +#ifndef OUTPUTS_RESTART_OPMD_HPP_ +#define OUTPUTS_RESTART_OPMD_HPP_ +//! \file restart_opmd.hpp +// \brief Provides support for restarting from OpenPMD output + +#include +#include +#include + +#include "basic_types.hpp" +#include "interface/swarm_default_names.hpp" +#include "openPMD/Iteration.hpp" +#include "openPMD/Series.hpp" +#include "outputs/restart.hpp" + +#include "mesh/domain.hpp" + +namespace parthenon { + +class Mesh; +class Param; + +class RestartReaderOPMD : public RestartReader { + public: + explicit RestartReaderOPMD(const char *filename); + + [[nodiscard]] SparseInfo GetSparseInfo() const override; + + [[nodiscard]] MeshInfo GetMeshInfo() const override; + + [[nodiscard]] SimTime GetTimeInfo() const override; + + [[nodiscard]] std::string GetInputString() const override { + return it->getAttribute("InputFile").get(); + }; + + // Return output format version number. Return -1 if not existent. + [[nodiscard]] int GetOutputFormatVersion() const override; + + // Current not supported + [[nodiscard]] int HasGhost() const override { return 0; }; + + public: + // Gets data for all blocks on current rank. + // Assumes blocks are contiguous + // fills internal data for given pointer + void ReadBlocks(const std::string &name, IndexRange range, + const OutputUtils::VarInfo &info, std::vector &dataVec, + int file_output_format_version, Mesh *pmesh) const override; + + // Gets the data from a swarm var on current rank. Assumes all + // blocks are contiguous. Fills dataVec based on shape from swarmvar + // metadata. + template + void ReadSwarmVar(const std::string &swarmname, const std::string &varname, + const std::size_t count, const std::size_t offset, const Metadata &m, + std::vector &data_vec) { + openPMD::ParticleSpecies swm = it->particles[swarmname]; + + const auto &shape = m.Shape(); + const int rank = shape.size(); + std::size_t nvar = 1; + for (int i = 0; i < rank; ++i) { + nvar *= shape[rank - 1 - i]; + } + std::size_t total_count = nvar * count; + if (data_vec.size() < total_count) { // greedy re-alloc + data_vec.resize(total_count); + } + + std::string particle_record; + std::string particle_record_component; + for (auto n = 0; n < nvar; n++) { + if (varname == swarm_position::x::name()) { + particle_record = "position"; + particle_record_component = "x"; + } else if (varname == swarm_position::y::name()) { + particle_record = "position"; + particle_record_component = "y"; + } else if (varname == swarm_position::z::name()) { + particle_record = "position"; + particle_record_component = "z"; + } else { + particle_record = varname; + particle_record_component = + rank == 0 ? openPMD::MeshRecordComponent::SCALAR : std::to_string(n); + } + + openPMD::RecordComponent rc = swm[particle_record][particle_record_component]; + rc.loadChunkRaw(&data_vec[n * count], {offset}, {count}); + } + + // Now actually read the registered chunks form disk + it->seriesFlush(); + } + + void ReadSwarmVar(const std::string &swarmname, const std::string &varname, + const std::size_t count, const std::size_t offset, const Metadata &m, + std::vector &dataVec) override { + ReadSwarmVar<>(swarmname, varname, count, offset, m, dataVec); + }; + void ReadSwarmVar(const std::string &swarmname, const std::string &varname, + const std::size_t count, const std::size_t offset, const Metadata &m, + std::vector &dataVec) override { + ReadSwarmVar<>(swarmname, varname, count, offset, m, dataVec); + }; + + // Gets the counts and offsets for MPI ranks for the meshblocks set + // by the indexrange. Returns the total count on this rank. + [[nodiscard]] std::size_t GetSwarmCounts(const std::string &swarm, + const IndexRange &range, + std::vector &counts, + std::vector &offsets) override; + + void ReadParams(const std::string &name, Params &p) override; + + template + void RestoreViewAttribute(const std::string &full_path, T &view) { + auto rank_and_dims = + it->getAttribute(full_path + ".rankdims").get>(); + // Resize view. + typename T::array_layout layout; + for (int d = 0; d < rank_and_dims[0]; ++d) { + layout.dimension[d] = rank_and_dims[1 + d]; + } + // Cannot use Kokkos::resize here as it's ambiguous at this point. + // Also, resize() interally also just create a new view. + view = T(Kokkos::view_alloc(Kokkos::WithoutInitializing, view.label()), layout); + auto view_h = Kokkos::create_mirror_view(HostMemSpace(), view); + + using base_t = typename std::remove_pointer::type; + auto flat_data = it->getAttribute(full_path).get>(); + for (auto i = 0; i < view_h.size(); i++) { + view_h.data()[i] = flat_data[i]; + } + Kokkos::deep_copy(view, view_h); + } + // closes out the restart file + // perhaps belongs in a destructor? + void Close(); + + private: + const std::string filename_; + + openPMD::Series series; + // Iteration is a pointer because it cannot be default constructed (it depends on the + // Series). + std::unique_ptr it; + + template + void ReadAllParamsOfType(const std::string &prefix, Params ¶ms); + template + void ReadAllParamsOfMultipleTypes(const std::string &prefix, Params &p); + template + void ReadAllParams(const std::string &pkg_name, Params &p); +}; + +} // namespace parthenon +#endif // OUTPUTS_RESTART_OPMD_HPP_ diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 26b5909e8775..ab2d3b29cf21 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -38,6 +38,7 @@ #include "outputs/output_utils.hpp" #include "outputs/restart.hpp" #include "outputs/restart_hdf5.hpp" +#include "outputs/restart_opmd.hpp" #include "utils/error_checking.hpp" #include "utils/utils.hpp" @@ -100,6 +101,8 @@ ParthenonStatus ParthenonManager::ParthenonInitEnv(int argc, char *argv[]) { // Read input from restart file if (fs::path(arg.restart_filename).extension() == ".rhdf") { restartReader = std::make_unique(arg.restart_filename); + } else if (fs::path(arg.restart_filename).extension() == ".bp") { + restartReader = std::make_unique(arg.restart_filename); } else { PARTHENON_FAIL("Unsupported restart file format."); } @@ -265,7 +268,8 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { // Currently supports versions 3 and 4. const auto file_output_format_ver = resfile.GetOutputFormatVersion(); - if (file_output_format_ver < HDF5::OUTPUT_VERSION_FORMAT - 1) { + // TODO(pgrete) figure out what to do about versions of different outputs + if (false && file_output_format_ver < HDF5::OUTPUT_VERSION_FORMAT - 1) { std::stringstream msg; msg << "File format version " << file_output_format_ver << " not supported. " << "Current format is " << HDF5::OUTPUT_VERSION_FORMAT << std::endl; @@ -325,7 +329,7 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { } // Read relevant data from the hdf file, this works for dense and sparse variables try { - resfile.ReadBlocks(label, myBlocks, v_info, tmp, file_output_format_ver); + resfile.ReadBlocks(label, myBlocks, v_info, tmp, file_output_format_ver, &rm); } catch (std::exception &ex) { std::cout << "[" << Globals::my_rank << "] WARNING: Failed to read variable " << label << " from restart file:" << std::endl @@ -355,7 +359,8 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { // Double note that this also needs to be update in case // we update the HDF5 infrastructure! - if (file_output_format_ver >= HDF5::OUTPUT_VERSION_FORMAT - 1) { + // TODO(pgrete) figure out what to do about versions of different outputs + if (true || file_output_format_ver >= HDF5::OUTPUT_VERSION_FORMAT - 1) { OutputUtils::PackOrUnpackVar( v_info, resfile.HasGhost() != 0, index, [&](auto index, int topo, int t, int u, int v, int k, int j, int i) { diff --git a/src/utils/mpi_types.hpp b/src/utils/mpi_types.hpp index a990cecc7e3c..6f5e117504ad 100644 --- a/src/utils/mpi_types.hpp +++ b/src/utils/mpi_types.hpp @@ -1,4 +1,8 @@ //======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2021-2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== // (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los @@ -34,6 +38,11 @@ inline MPI_Datatype MPITypeMap::type() { return MPI_PARTHENON_REAL; } +template <> +inline MPI_Datatype MPITypeMap::type() { + return MPI_INT64_T; +} + template <> inline MPI_Datatype MPITypeMap::type() { return MPI_INT; @@ -44,6 +53,12 @@ inline MPI_Datatype MPITypeMap::type() { return MPI_CXX_BOOL; } +template <> +inline MPI_Datatype MPITypeMap::type() { + // TODO(pgrete) do we need special checks here wrt to conflicts on MacOS? + return MPI_UINT64_T; +} + } // namespace parthenon #endif diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt index d63f81a03ffc..79ddd662f3e7 100644 --- a/tst/regression/CMakeLists.txt +++ b/tst/regression/CMakeLists.txt @@ -152,6 +152,21 @@ if (ENABLE_HDF5) endif() +if (PARTHENON_ENABLE_OPENPMD) + + # h5py is needed for restart and hdf5 test + list(APPEND REQUIRED_PYTHON_MODULES openpmd_api) + + # Restart + list(APPEND TEST_DIRS restart_opmd) + list(APPEND TEST_PROCS ${NUM_MPI_PROC_TESTING}) + list(APPEND TEST_ARGS "--driver ${PROJECT_BINARY_DIR}/example/advection/advection-example \ + --driver_input ${CMAKE_CURRENT_SOURCE_DIR}/test_suites/restart_opmd/parthinput.restart \ + --num_steps 2") + list(APPEND EXTRA_TEST_LABELS "") + + endif() + # Any external modules that are required by python can be added to REQUIRED_PYTHON_MODULES # list variable, before including TestSetup.cmake. list(APPEND REQUIRED_PYTHON_MODULES numpy) diff --git a/tst/regression/test_suites/restart_opmd/__init__.py b/tst/regression/test_suites/restart_opmd/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/tst/regression/test_suites/restart_opmd/parthinput.restart b/tst/regression/test_suites/restart_opmd/parthinput.restart new file mode 100644 index 000000000000..768453160a7b --- /dev/null +++ b/tst/regression/test_suites/restart_opmd/parthinput.restart @@ -0,0 +1,60 @@ +# ======================================================================================== +# Parthenon performance portable AMR framework +# Copyright(C) 2024 The Parthenon collaboration +# Licensed under the 3-clause BSD License, see LICENSE file for details +# ======================================================================================== + + +problem_id = advection + + +refinement = adaptive +numlevel = 3 + +nx1 = 64 +x1min = -0.5 +x1max = 0.5 +ix1_bc = periodic +ox1_bc = periodic + +nx2 = 64 +x2min = -0.5 +x2max = 0.5 +ix2_bc = periodic +ox2_bc = periodic + +nx3 = 1 +x3min = -0.5 +x3max = 0.5 +ix3_bc = periodic +ox3_bc = periodic + + +nx1 = 16 +nx2 = 16 +nx3 = 1 + + +nlim = -1 +tlim = 0.2 +integrator = rk2 +ncycle_out_mesh = -10000 + + +cfl = 0.45 +vx = 1.0 +vy = 1.0 +vz = 1.0 +profile = hard_sphere + +refine_tol = 0.3 # control the package specific refinement tagging function +derefine_tol = 0.03 +compute_error = false +num_vars = 1 # number of variables +vec_size = 1 # size of each variable +fill_derived = false # whether to fill one-copy test vars + + +file_type = openpmd +dt = 0.050 + diff --git a/tst/regression/test_suites/restart_opmd/parthinput_override.restart b/tst/regression/test_suites/restart_opmd/parthinput_override.restart new file mode 100644 index 000000000000..c5b368aebcbd --- /dev/null +++ b/tst/regression/test_suites/restart_opmd/parthinput_override.restart @@ -0,0 +1,9 @@ +# ======================================================================================== +# Parthenon performance portable AMR framework +# Copyright(C) 2024 The Parthenon collaboration +# Licensed under the 3-clause BSD License, see LICENSE file for details +# ======================================================================================== + +# Testing to override parameters in a restart file from an input file + +problem_id=silver diff --git a/tst/regression/test_suites/restart_opmd/restart_opmd.py b/tst/regression/test_suites/restart_opmd/restart_opmd.py new file mode 100644 index 000000000000..3548418f61d4 --- /dev/null +++ b/tst/regression/test_suites/restart_opmd/restart_opmd.py @@ -0,0 +1,161 @@ +# ======================================================================================== +# Parthenon performance portable AMR framework +# Copyright(C) 2024 The Parthenon collaboration +# Licensed under the 3-clause BSD License, see LICENSE file for details +# ======================================================================================== + +# Modules +import sys +import utils.test_case +import numpy as np + + +# To prevent littering up imported folders with .pyc files or __pycache_ folder +sys.dont_write_bytecode = True + + +class TestCase(utils.test_case.TestCaseAbs): + def Prepare(self, parameters, step): + # enable coverage testing on pass where restart + # files are both read and written + parameters.coverage_status = "both" + + # run baseline (to the very end) + if step == 1: + parameters.driver_cmd_line_args = ["parthenon/job/problem_id=gold"] + # restart from an early snapshot + elif step == 2: + # TODO(pgrete or someone else) ideally we want to restart from a later snapshot + # BUT results are not bitwise identical for AMR runs. PG thinks this is + # related to not storing the deref counter (and similar) and also thinks + # it's worth fixing. + parameters.driver_cmd_line_args = [ + "-r", + "gold.out1.00000.bp", + "-i", + f"{parameters.parthenon_path}/tst/regression/test_suites/restart_opmd/parthinput_override.restart", + ] + + return parameters + + def Analyse(self, parameters): + try: + import openpmd_api as opmd + except ModuleNotFoundError: + print("Couldn't find required openpmd_api module to compare test results.") + return False + success = True + + def compare_attributes(series_a, series_b): + skip_attributes = [ + "iterationFormat", # Stores the file name format. Expected to differ. + "WallTime", + "InputFile", # Is updated during runtime, e.g., startime and thus differs + ] + all_equal = True + for attr in series_a.attributes: + if series_b.contains_attribute(attr): + attr_a = series_a.get_attribute(attr) + attr_b = series_b.get_attribute(attr) + if attr not in skip_attributes and attr_a != attr_b: + print( + f"Mismatch in attribute '{attr}'. " + f"'{attr_a}' versus '{attr_b}'\n" + ) + all_equal = False + else: + print(f"Missing attribute '{attr}' in second file.") + all_equal = False + return all_equal + + # need series in order to flush + def compare_data(it_a, it_b, series_a, series_b): + all_equal = True + for mesh_name, mesh_a in it_a.meshes.items(): + if mesh_name not in it_b.meshes: + print(f"Missing mesh '{mesh_name}' in second file.") + all_equal = False + continue + mesh_b = it_b.meshes[mesh_name] + + for comp_name, comp_a in mesh_a.items(): + if comp_name not in mesh_b: + print( + f"Missing component '{comp_name}' in mesh '{mesh_name}' of second file." + ) + all_equal = False + continue + comp_b = mesh_b[comp_name] + + if comp_a.shape != comp_b.shape: + print( + f"Mismatch is mech record component shapes of " + " compontent '{comp_name}' in mesh '{mesh_name}': " + f"{comp_a.shape} versus {comp_b.shape}\n" + ) + all_equal = False + continue + + # Given that the shapes are guaranteed to match (follow the check above) + # we can load chunks from both files. + # Note that we have to go over chunks as data might be sparse on disk so + # loading the entire record will contain gargabe in sparse places. + data_a = np.empty(comp_a.shape) + data_a[:] = np.nan + data_b = np.copy(data_a) + for chunk in comp_a.available_chunks(): + # Following OpenPMD-viewer `chunk_to_slice` here + # https://github.com/openPMD/openPMD-viewer/blob/6eccb608893d2c9b8d158d950c3f0451898a80f6/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/utilities.py#L14 + stops = [a + b for a, b in zip(chunk.offset, chunk.extent)] + indices_per_dim = zip(chunk.offset, stops) + sl = tuple( + map(lambda s: slice(s[0], s[1], None), indices_per_dim) + ) + + tmp = comp_a[sl] + series_a.flush() + data_a[sl] = tmp + + tmp = comp_b[sl] + series_b.flush() + data_b[sl] = tmp + + try: + np.testing.assert_array_max_ulp(data_a, data_b) + except AssertionError as err: + print( + f"Data of component '{comp_name}' in mesh '{mesh_name}' does not match:\n" + f"{err}\n" + ) + all_equal = False + continue + + return all_equal + + def compare_files(idx_it): + all_good = True + series_gold = opmd.Series("gold.out1.%T.bp/", opmd.Access.read_only) + series_silver = opmd.Series("silver.out1.%T.bp/", opmd.Access.read_only) + + # PG: yes, this is inefficient but keeps the logic simple + all_good &= compare_attributes(series_gold, series_silver) + all_good &= compare_attributes(series_silver, series_gold) + + it_gold = series_gold.iterations[idx_it] + it_silver = series_silver.iterations[idx_it] + all_good &= compare_attributes(it_gold, it_silver) + all_good &= compare_attributes(it_silver, it_gold) + + all_good &= compare_data(it_silver, it_gold, series_silver, series_gold) + all_good &= compare_data(it_gold, it_silver, series_gold, series_silver) + + return all_good + + # comapre a few files throughout the simulations + success &= compare_files(1) + success &= compare_files(2) + success &= compare_files(3) + success &= compare_files(4) + # success &= compare_files("final") + + return success diff --git a/tst/unit/test_unit_params.cpp b/tst/unit/test_unit_params.cpp index 633d7f856163..aae19bd1916c 100644 --- a/tst/unit/test_unit_params.cpp +++ b/tst/unit/test_unit_params.cpp @@ -1,6 +1,6 @@ //======================================================================================== -// Athena++ astrophysical MHD code -// Copyright(C) 2014 James M. Stone and other code contributors +// Parthenon performance portable AMR framework +// Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== // (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved. @@ -16,6 +16,8 @@ //======================================================================================== #include +#include +#include #include #include @@ -24,7 +26,12 @@ #include "config.hpp" #include "interface/params.hpp" #include "kokkos_abstraction.hpp" +#include "openPMD/Series.hpp" #include "outputs/parthenon_hdf5.hpp" +#include "outputs/parthenon_opmd.hpp" +#include "outputs/restart_hdf5.hpp" +#include "outputs/restart_opmd.hpp" +#include "parthenon_array_generic.hpp" using parthenon::Params; using parthenon::Real; @@ -136,9 +143,22 @@ TEST_CASE("when hasKey is called", "[hasKey]") { } } -#ifdef ENABLE_HDF5 - -TEST_CASE("A set of params can be dumped to file", "[params][output]") { +#if defined(ENABLE_HDF5) && defined(PARTHENON_ENABLE_OPENPMD) +using parthenon::RestartReaderHDF5; +using parthenon::RestartReaderOPMD; +using OutputTypes = std::tuple; +#elif defined(ENABLE_HDF5) +using parthenon::RestartReaderHDF5; +using OutputTypes = std::tuple; +#elif defined(PARTHENON_ENABLE_OPENPMD) +using parthenon::RestartReaderOPMD; +using OutputTypes = std::tuple; +#else +using OutputTypes = std::tuple<>; +#endif + +TEMPLATE_LIST_TEST_CASE("A set of params can be dumped to file", "[params][output]", + OutputTypes) { GIVEN("A params object with a few kinds of objects") { Params params; const auto restart = Params::Mutability::Restart; @@ -163,43 +183,97 @@ TEST_CASE("A set of params can be dumped to file", "[params][output]") { Kokkos::deep_copy(arr2d, arr2d_h); params.Add("arr2d", arr2d); - parthenon::HostArray2D hostarr("hostarr2d", 2, 3); + // "Vectors" of bools have some special sauce under the hood so let's try the logic + // with a plain view + Kokkos::View bool1d("boolview", 10); + auto bool1d_h = Kokkos::create_mirror_view(bool1d); + for (int i = 0; i < 10; ++i) { + bool1d_h(i) = i % 2; + } + Kokkos::deep_copy(bool1d, bool1d_h); + params.Add("bool1d", bool1d); + + parthenon::HostArray2D hostarr2d("hostarr2d", 2, 3); for (int i = 0; i < 2; ++i) { for (int j = 0; j < 3; ++j) { - hostarr(i, j) = 2 * i + j + 1; + hostarr2d(i, j) = 2 * i + j + 1; } } - params.Add("hostarr2d", hostarr, restart); + params.Add("hostarr2d", hostarr2d, restart); - THEN("We can output to hdf5") { - const std::string filename = "params_test.h5"; - const std::string groupname = "params"; + THEN("We can output") { + std::string filename; + const std::string groupname = "Params"; const std::string prefix = "test_pkg"; - using namespace parthenon::HDF5; - { + if constexpr (std::is_same_v) { + using namespace parthenon::HDF5; + filename = "params_test.h5"; + H5F file = H5F::FromHIDCheck( H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)); auto group = MakeGroup(file, groupname); params.WriteAllToHDF5(prefix, group); + } else if constexpr (std::is_same_v) { + filename = ("params_test.%05T.bp"); + auto series = openPMD::Series(filename, openPMD::Access::CREATE); + series.setIterationEncoding(openPMD::IterationEncoding::fileBased); + auto it = series.iterations[0]; + parthenon::OpenPMDUtils::WriteAllParams(params, prefix, &it); + } else { + FAIL("This logic is flawed. I should not be here."); } - AND_THEN("We can directly read the relevant data from the hdf5 file") { - H5F file = - H5F::FromHIDCheck(H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); - const H5O obj = H5O::FromHIDCheck(H5Oopen(file, groupname.c_str(), H5P_DEFAULT)); - + AND_THEN("We can directly read the relevant data from the file") { Real in_scalar; - HDF5ReadAttribute(obj, prefix + "/scalar", in_scalar); - REQUIRE(std::abs(scalar - in_scalar) <= 1e-10); - std::vector in_vector; - HDF5ReadAttribute(obj, prefix + "/vector", in_vector); + // deliberately the wrong size + parthenon::ParArray2D in_arr2d("myarr", 1, 1); + parthenon::HostArray2D in_hostarr2d("hostarr2d", 2, 3); + Kokkos::View in_bool1d("in_bool1d", 5); + + if constexpr (std::is_same_v) { + H5F file = + H5F::FromHIDCheck(H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); + const H5O obj = + H5O::FromHIDCheck(H5Oopen(file, groupname.c_str(), H5P_DEFAULT)); + + HDF5ReadAttribute(obj, prefix + "/scalar", in_scalar); + HDF5ReadAttribute(obj, prefix + "/vector", in_vector); + HDF5ReadAttribute(obj, prefix + "/arr2d", in_arr2d); + HDF5ReadAttribute(obj, prefix + "/hostarr2d", in_hostarr2d); + HDF5ReadAttribute(obj, prefix + "/bool1d", in_bool1d); + } else if constexpr (std::is_same_v) { + auto series = openPMD::Series(filename, openPMD::Access::READ_ONLY); + auto it = series.iterations[0]; + // Note that we're explicitly using `delim` here which tests the character + // replacement of '/' in the WriteAllParams function. + using parthenon::OpenPMDUtils::delim; + + in_scalar = + it.getAttribute(groupname + delim + prefix + delim + "scalar").get(); + + in_vector = it.getAttribute(groupname + delim + prefix + delim + "vector") + .get>(); + + // Technically, we're not reading "directly" here but the restart reader ctor + // literally just opens the file. + auto resfile = RestartReaderOPMD(filename.c_str()); + auto &in_arr2d_v = in_arr2d.KokkosView(); + resfile.RestoreViewAttribute(groupname + delim + prefix + delim + "arr2d", + in_arr2d_v); + + auto &in_hostarr2d_v = in_hostarr2d.KokkosView(); + resfile.RestoreViewAttribute(groupname + delim + prefix + delim + "hostarr2d", + in_hostarr2d_v); + // TODO(pgrete) make this work and also add checks for correctness below + // resfile.RestoreViewAttribute(groupname + delim + prefix + delim + "bool1d", + // in_bool1d); + } + REQUIRE(scalar == in_scalar); + for (int i = 0; i < vector.size(); ++i) { REQUIRE(in_vector[i] == vector[i]); } - // deliberately the wrong size - parthenon::ParArray2D in_arr2d("myarr", 1, 1); - HDF5ReadAttribute(obj, prefix + "/arr2d", in_arr2d); REQUIRE(in_arr2d.extent_int(0) == arr2d.extent_int(0)); REQUIRE(in_arr2d.extent_int(1) == arr2d.extent_int(1)); int nwrong = 1; @@ -211,9 +285,17 @@ TEST_CASE("A set of params can be dumped to file", "[params][output]") { }, nwrong); REQUIRE(nwrong == 0); + + REQUIRE(in_hostarr2d.extent_int(0) == hostarr2d.extent_int(0)); + REQUIRE(in_hostarr2d.extent_int(1) == hostarr2d.extent_int(1)); + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 3; ++j) { + REQUIRE(hostarr2d(i, j) == in_hostarr2d(i, j)); + } + } } - AND_THEN("We can restart a params object from the HDF5 file") { + AND_THEN("We can restart a params object from the file") { Params rparams; // init the params object to restart into @@ -232,24 +314,30 @@ TEST_CASE("A set of params can be dumped to file", "[params][output]") { parthenon::HostArray2D test_hostarr("hostarr2d", 1, 1); rparams.Add("hostarr2d", test_hostarr, restart); - H5F file = - H5F::FromHIDCheck(H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); - const H5G obj = H5G::FromHIDCheck(H5Oopen(file, groupname.c_str(), H5P_DEFAULT)); - rparams.ReadFromRestart(prefix, obj); + if constexpr (std::is_same_v) { + H5F file = + H5F::FromHIDCheck(H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); + const H5G obj = + H5G::FromHIDCheck(H5Oopen(file, groupname.c_str(), H5P_DEFAULT)); + rparams.ReadFromRestart(prefix, obj); + } else if constexpr (std::is_same_v) { + auto resfile = RestartReaderOPMD(filename.c_str()); + resfile.ReadParams(prefix, rparams); + } AND_THEN("The values for the restartable params are updated to match the file") { auto test_scalar = rparams.Get("scalar"); - REQUIRE(std::abs(test_scalar - scalar) <= 1e-10); + REQUIRE(test_scalar == scalar); auto test_bool = rparams.Get("boolscalar"); REQUIRE(test_bool == boolscalar); - auto test_hostarr = params.Get>("hostarr2d"); - REQUIRE(test_hostarr.extent_int(0) == hostarr.extent_int(0)); - REQUIRE(test_hostarr.extent_int(1) == hostarr.extent_int(1)); - for (int i = 0; i < hostarr.extent_int(0); ++i) { - for (int j = 0; j < hostarr.extent_int(1); ++j) { - REQUIRE(test_hostarr(i, j) == hostarr(i, j)); + auto test_hostarr = rparams.Get>("hostarr2d"); + REQUIRE(test_hostarr.extent_int(0) == hostarr2d.extent_int(0)); + REQUIRE(test_hostarr.extent_int(1) == hostarr2d.extent_int(1)); + for (int i = 0; i < hostarr2d.extent_int(0); ++i) { + for (int j = 0; j < hostarr2d.extent_int(1); ++j) { + REQUIRE(test_hostarr(i, j) == hostarr2d(i, j)); } } } @@ -264,5 +352,3 @@ TEST_CASE("A set of params can be dumped to file", "[params][output]") { } } } - -#endif // ENABLE_HDF5