diff --git a/docs/source/usage/plugins/phaseSpace.rst b/docs/source/usage/plugins/phaseSpace.rst index 63a73c2a32a..e3fc40f7250 100644 --- a/docs/source/usage/plugins/phaseSpace.rst +++ b/docs/source/usage/plugins/phaseSpace.rst @@ -8,7 +8,7 @@ This plugin creates a 2D phase space image for a user-given spatial and momentum External Dependencies ^^^^^^^^^^^^^^^^^^^^^ -The plugin is available as soon as the :ref:`libSplash and HDF5 libraries ` are compiled in. +The plugin is available as soon as the :ref:`openPMD API ` is compiled in. .cfg file ^^^^^^^^^ @@ -19,13 +19,13 @@ Example for *y-pz* phase space for the *electron* species (``.cfg`` file macro): # Calculate a 2D phase space # - momentum range in m_e c - TGB_ePSypz="--e_phaseSpace.period 10 --e_phaseSpace.filter all --e_phaseSpace.space y --e_phaseSpace.momentum pz --e_phaseSpace.min -1.0 --e_phaseSpace.max 1.0" + TGB_ePSypz="--e_phaseSpace.period 10 --e_phaseSpace.filter all --e_phaseSpace.space y --e_phaseSpace.momentum pz --e_phaseSpace.min -1.0 --e_phaseSpace.max 1.0 --e_phaseSpace.ext h5" The distinct options are (assuming a species ``e`` for electrons): ====================================== ======================================================== ============================ -Option Usage Unit +Option Usage Unit ====================================== ======================================================== ============================ ``--e_phaseSpace.period `` calculate each N steps *none* ``--e_phaseSpace.filter`` Use filtered particles. Available filters are set up in *none* @@ -34,6 +34,7 @@ Option Usage ``--e_phaseSpace.momentum `` momentum coordinate of the 2D phase space *none* ``--e_phaseSpace.min `` minimum of the momentum range :math:`m_\mathrm{species} c` ``--e_phaseSpace.max `` maximum of the momentum range :math:`m_\mathrm{species} c` +``--e_phaseSpace.ext `` filename extension for openPMD backend *none* ====================================== ======================================================== ============================ Memory Complexity @@ -52,11 +53,23 @@ negligible. Output ^^^^^^ -The 2D histograms are stored in ``.hdf5`` files in the ``simOutput/phaseSpace/`` directory. +The 2D histograms are stored in the ``simOutput/phaseSpace/`` directory, by default in ``.h5`` files. A file is created per species, phasespace selection and time step. Values are given as *charge density* per phase space bin. -In order to scale to a simpler *charge of particles* per :math:`\mathrm{d}r_i` and :math:`\mathrm{d}p_i` -bin multiply by the cell volume ``dV``. +In order to scale to a simpler *charge of particles* per :math:`\mathrm{d}r_i` and :math:`\mathrm{d}p_i` -bin multiply by the cell volume ``dV`` (written as an attribute of the openPMD Mesh). + +The output writes a number of non-standard custom openPMD attributes: + +* ``p_min`` and ``p_max``: The lower and upper bounds for the momentum axis, respectively. +* ``dr``: The spacing of the spatial axis in PIConGPU units. +* ``dV``: The volume of a phase space cell. Relates to ``dr`` via ``dV = dp * dr`` where ``dp`` would be the grid spacing along the momentum axis. +* ``dr_unit``: The SI scaling for the spatial axis. Use this instead of ``gridUnitSI``. +* ``p_unit``: The SI scaling for the momentum axis. Use this instead of ``gridUnitSI``. +* ``globalDomainOffset``, ``globalDomainSize`` and ``globalDomainAxisLabels``: Information on the global domain. +* ``totalDomainOffset``, ``totalDomainSize`` and ``totalDomainAxisLabels``: Information on the total domain. + Please consult the `PIConGPU wiki `_ for explanations on the meaning of global and total domain. +* ``sim_unit``: SI scaling for the charge density values. Alias for ``unitSI``. Analysis Tools ^^^^^^^^^^^^^^ @@ -223,7 +236,8 @@ Known Limitations - charge deposition uses the counter shape for now (would need one more write to neighbors to evaluate it correctly according to the shape) - the user has to define the momentum range in advance - the resolution is fixed to ``1024 bins`` in momentum and the number of cells in the selected spatial dimension -- this plugin does not yet use :ref:`openPMD markup `. +- While the openPMD standard `has already been updated `_ to support phase space data, the openPMD API does not yet implement this part. + The openPMD attribute ``gridUnitSI`` and ``gridUnitDimension`` can hence not be correctly written yet and should be ignored in favor of the custom attributes written by this plugin. References ^^^^^^^^^^ diff --git a/include/picongpu/plugins/PhaseSpace/AxisDescription.hpp b/include/picongpu/plugins/PhaseSpace/AxisDescription.hpp index 71f67dab0eb..1fee9621092 100644 --- a/include/picongpu/plugins/PhaseSpace/AxisDescription.hpp +++ b/include/picongpu/plugins/PhaseSpace/AxisDescription.hpp @@ -19,6 +19,8 @@ #pragma once +#include + namespace picongpu { /** 2D Phase Space Selection @@ -47,6 +49,36 @@ namespace picongpu y = 1u, z = 2u }; + + std::string momentumAsString() const + { + switch(momentum) + { + case px: + return "px"; + case py: + return "py"; + case pz: + return "pz"; + default: + throw std::runtime_error("Unreachable!"); + } + } + + std::string spaceAsString() const + { + switch(space) + { + case x: + return "x"; + case y: + return "y"; + case z: + return "z"; + default: + throw std::runtime_error("Unreachable!"); + } + } }; } /* namespace picongpu */ diff --git a/include/picongpu/plugins/PhaseSpace/DumpHBufferOpenPMD.hpp b/include/picongpu/plugins/PhaseSpace/DumpHBufferOpenPMD.hpp new file mode 100644 index 00000000000..8de5ab8ca9b --- /dev/null +++ b/include/picongpu/plugins/PhaseSpace/DumpHBufferOpenPMD.hpp @@ -0,0 +1,224 @@ +/* Copyright 2013-2021 Axel Huebl, Rene Widera + * + * This file is part of PIConGPU. + * + * PIConGPU is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PIConGPU is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PIConGPU. + * If not, see . + */ + +#pragma once + +#include "picongpu/simulation_defines.hpp" + +#include "picongpu/plugins/PhaseSpace/AxisDescription.hpp" +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace picongpu +{ + class DumpHBuffer + { + private: + using SuperCellSize = typename MappingDesc::SuperCellSize; + + public: + /** Dump the PhaseSpace host Buffer + * + * \tparam Type the HBuffers element type + * \tparam int the HBuffers dimension + * \param hBuffer const reference to the hBuffer, including guard cells in spatial dimension + * \param axis_element plot to create: e.g. py, x from momentum/spatial-coordinate + * \param unit sim unit of the buffer + * \param strSpecies unique short hand name of the species + * \param filenameSuffix infix + extension part of openPMD filename + * \param currentStep current time step + * \param mpiComm communicator of the participating ranks + */ + template + void operator()( + const pmacc::container::HostBuffer& hBuffer, + const AxisDescription axis_element, + const std::pair axis_p_range, + const float_64 pRange_unit, + const float_64 unit, + const std::string strSpecies, + const std::string filenameExtension, + const std::string jsonConfig, + const uint32_t currentStep, + MPI_Comm mpiComm) const + { + using Type = T_Type; + + /** file name ***************************************************** + * phaseSpace/PhaseSpace_xpy_timestep.h5 */ + std::string fCoords("xyz"); + std::ostringstream openPMDFilename; + openPMDFilename << "phaseSpace/PhaseSpace_" << strSpecies << "_" << fCoords.at(axis_element.space) << "p" + << fCoords.at(axis_element.momentum) << "_%T." << filenameExtension; + + /** get size of the fileWriter communicator ***********************/ + int size; + MPI_CHECK(MPI_Comm_size(mpiComm, &size)); + + /** create parallel domain collector ******************************/ + ::openPMD::Series series(openPMDFilename.str(), ::openPMD::Access::CREATE, jsonConfig); + ::openPMD::Iteration iteration = series.iterations[currentStep]; + + const std::string software("PIConGPU"); + + std::stringstream softwareVersion; + softwareVersion << PICONGPU_VERSION_MAJOR << "." << PICONGPU_VERSION_MINOR << "." + << PICONGPU_VERSION_PATCH; + if(!std::string(PICONGPU_VERSION_LABEL).empty()) + softwareVersion << "-" << PICONGPU_VERSION_LABEL; + series.setSoftware(software, softwareVersion.str()); + + pmacc::GridController& gc = pmacc::Environment::get().GridController(); + + /** calculate GUARD offset in the source hBuffer *****************/ + const uint32_t rGuardCells + = SuperCellSize().toRT()[axis_element.space] * GuardSize::toRT()[axis_element.space]; + + /** calculate local and global size of the phase space ***********/ + const uint32_t numSlides = MovingWindow::getInstance().getSlideCounter(currentStep); + const SubGrid& subGrid = Environment::get().SubGrid(); + const std::uint64_t rLocalOffset = subGrid.getLocalDomain().offset[axis_element.space]; + const std::uint64_t rLocalSize = int(hBuffer.size().y() - 2 * rGuardCells); + const std::uint64_t rGlobalSize = subGrid.getGlobalDomain().size[axis_element.space]; + PMACC_VERIFY(int(rLocalSize) == subGrid.getLocalDomain().size[axis_element.space]); + + /* globalDomain of the phase space */ + ::openPMD::Extent globalPhaseSpace_extent{rGlobalSize, hBuffer.size().x()}; + + /* global moving window meta information */ + ::openPMD::Offset globalPhaseSpace_offset{0, 0}; + std::uint64_t globalMovingWindowOffset = 0; + std::uint64_t globalMovingWindowSize = rGlobalSize; + if(axis_element.space == AxisDescription::y) /* spatial axis == y */ + { + globalPhaseSpace_offset[0] = numSlides * rLocalSize; + Window window = MovingWindow::getInstance().getWindow(currentStep); + globalMovingWindowOffset = window.globalDimensions.offset[axis_element.space]; + globalMovingWindowSize = window.globalDimensions.size[axis_element.space]; + } + + /* localDomain: offset of it in the globalDomain and size */ + ::openPMD::Offset localPhaseSpace_offset{rLocalOffset, 0}; + ::openPMD::Extent localPhaseSpace_extent{rLocalSize, hBuffer.size().x()}; + + /** Dataset Name **************************************************/ + std::ostringstream dataSetName; + /* xpx or ypz or ... */ + dataSetName << strSpecies << "_" << fCoords.at(axis_element.space) << "p" + << fCoords.at(axis_element.momentum); + + /** debug log *****************************************************/ + int rank; + MPI_CHECK(MPI_Comm_rank(mpiComm, &rank)); + { + std::stringstream offsetAsString, localExtentAsString, globalExtentAsString; + offsetAsString << "[" << localPhaseSpace_offset[0] << ", " << localPhaseSpace_offset[1] << "]"; + localExtentAsString << "[" << localPhaseSpace_extent[0] << ", " << localPhaseSpace_extent[1] << "]"; + globalExtentAsString << "[" << globalPhaseSpace_extent[0] << ", " << globalPhaseSpace_extent[1] << "]"; + log( + "Dump buffer %1% to %2% at offset %3% with size %4% for total size %5% for rank %6% / %7%") + % (*(hBuffer.origin()(0, rGuardCells))) % dataSetName.str() % offsetAsString.str() + % localExtentAsString.str() % globalExtentAsString.str() % rank % size; + } + + /** write local domain ********************************************/ + + // avoid deadlock between not finished pmacc tasks and mpi calls in HDF5 + __getTransactionEvent().waitForFinished(); + + ::openPMD::Mesh mesh = iteration.meshes[dataSetName.str()]; + ::openPMD::MeshRecordComponent dataset = mesh[::openPMD::RecordComponent::SCALAR]; + + dataset.resetDataset({::openPMD::determineDatatype(), globalPhaseSpace_extent}); + std::shared_ptr data(&(*hBuffer.origin()(0, rGuardCells)), [](auto const&) {}); + dataset.storeChunk(data, localPhaseSpace_offset, localPhaseSpace_extent); + + /** meta attributes for the data set: unit, range, moving window **/ + + pmacc::Selection globalDomain = subGrid.getGlobalDomain(); + pmacc::Selection totalDomain = subGrid.getTotalDomain(); + // convert things to std::vector<> for the openPMD API to enjoy + std::vector globalDomainSize{&globalDomain.size[0], &globalDomain.size[0] + simDim}; + std::vector globalDomainOffset{&globalDomain.offset[0], &globalDomain.offset[0] + simDim}; + std::vector totalDomainSize{&totalDomain.size[0], &totalDomain.size[0] + simDim}; + std::vector totalDomainOffset{&totalDomain.offset[0], &totalDomain.offset[0] + simDim}; + std::vector globalDomainAxisLabels; + if(simDim == DIM2) + { + globalDomainAxisLabels = {"y", "x"}; // 2D: F[y][x] + } + if(simDim == DIM3) + { + globalDomainAxisLabels = {"z", "y", "x"}; // 3D: F[z][y][x] + } + + float_X const dr = cellSize[axis_element.space]; + + mesh.setAttribute("globalDomainSize", globalDomainSize); + mesh.setAttribute("globalDomainOffset", globalDomainOffset); + mesh.setAttribute("totalDomainSize", totalDomainSize); + mesh.setAttribute("totalDomainOffset", totalDomainOffset); + mesh.setAttribute("globalDomainAxisLabels", globalDomainAxisLabels); + mesh.setAttribute("totalDomainAxisLabels", globalDomainAxisLabels); + mesh.setAttribute("_global_start", globalPhaseSpace_offset); + mesh.setAttribute("_global_size", globalPhaseSpace_extent); + mesh.setAxisLabels({axis_element.spaceAsString(), axis_element.momentumAsString()}); + mesh.setAttribute("sim_unit", unit); + dataset.setUnitSI(unit); + { + using UD = ::openPMD::UnitDimension; + mesh.setUnitDimension({{UD::I, 1.0}, {UD::T, 1.0}, {UD::L, -1.0}}); // charge density + } + mesh.setAttribute("p_unit", pRange_unit); + mesh.setAttribute("p_min", axis_p_range.first); + mesh.setAttribute("p_max", axis_p_range.second); + mesh.setGridGlobalOffset({globalMovingWindowOffset * dr, axis_p_range.first}); + mesh.setAttribute("movingWindowOffset", globalMovingWindowOffset); + mesh.setAttribute("movingWindowSize", globalMovingWindowSize); + mesh.setAttribute("dr", dr); + mesh.setAttribute("dV", CELL_VOLUME); + mesh.setGridSpacing(std::vector{dr, CELL_VOLUME / dr}); + mesh.setAttribute("dr_unit", UNIT_LENGTH); + iteration.setDt(DELTA_T); + iteration.setTimeUnitSI(UNIT_TIME); + /* + * The value represents an aggregation over one cell, so any value is correct for the mesh position. + * Just use the center. + */ + dataset.setPosition(std::vector{0.5, 0.5}); + + /** close file ****************************************************/ + iteration.close(); + } + }; + +} /* namespace picongpu */ diff --git a/include/picongpu/plugins/PhaseSpace/DumpHBufferSplashP.hpp b/include/picongpu/plugins/PhaseSpace/DumpHBufferSplashP.hpp deleted file mode 100644 index 19562e043cf..00000000000 --- a/include/picongpu/plugins/PhaseSpace/DumpHBufferSplashP.hpp +++ /dev/null @@ -1,209 +0,0 @@ -/* Copyright 2013-2021 Axel Huebl, Rene Widera - * - * This file is part of PIConGPU. - * - * PIConGPU is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * PIConGPU is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with PIConGPU. - * If not, see . - */ - -#pragma once - -#include "picongpu/simulation_defines.hpp" - -#include "picongpu/traits/SplashToPIC.hpp" -#include "picongpu/traits/PICToSplash.hpp" - -#include "picongpu/plugins/PhaseSpace/AxisDescription.hpp" -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace picongpu -{ - class DumpHBuffer - { - private: - typedef typename MappingDesc::SuperCellSize SuperCellSize; - - public: - /** Dump the PhaseSpace host Buffer - * - * \tparam Type the HBuffers element type - * \tparam int the HBuffers dimension - * \param hBuffer const reference to the hBuffer, including guard cells in spatial dimension - * \param axis_element plot to create: e.g. py, x from momentum/spatial-coordinate - * \param unit sim unit of the buffer - * \param strSpecies unique short hand name of the species - * \param currentStep current time step - * \param mpiComm communicator of the participating ranks - */ - template - void operator()( - const pmacc::container::HostBuffer& hBuffer, - const AxisDescription axis_element, - const std::pair axis_p_range, - const float_64 pRange_unit, - const float_64 unit, - const std::string strSpecies, - const uint32_t currentStep, - MPI_Comm mpiComm) const - { - using namespace splash; - typedef T_Type Type; - const int bufDim = T_bufDim; - - /** file name ***************************************************** - * phaseSpace/PhaseSpace_xpy_timestep.h5 */ - std::string fCoords("xyz"); - std::ostringstream filename; - filename << "phaseSpace/PhaseSpace_" << strSpecies << "_" << fCoords.at(axis_element.space) << "p" - << fCoords.at(axis_element.momentum); - - /** get size of the fileWriter communicator ***********************/ - int size; - MPI_CHECK(MPI_Comm_size(mpiComm, &size)); - - /** create parallel domain collector ******************************/ - ParallelDomainCollector pdc(mpiComm, MPI_INFO_NULL, Dimensions(size, 1, 1), 10); - - pmacc::GridController& gc = pmacc::Environment::get().GridController(); - DataCollector::FileCreationAttr fAttr; - Dimensions mpiPosition(gc.getPosition()[axis_element.space], 0, 0); - fAttr.mpiPosition.set(mpiPosition); - - DataCollector::initFileCreationAttr(fAttr); - - pdc.open(filename.str().c_str(), fAttr); - - /** calculate GUARD offset in the source hBuffer *****************/ - const uint32_t rGuardCells - = SuperCellSize().toRT()[axis_element.space] * GuardSize::toRT()[axis_element.space]; - - /** calculate local and global size of the phase space ***********/ - const uint32_t numSlides = MovingWindow::getInstance().getSlideCounter(currentStep); - const SubGrid& subGrid = Environment::get().SubGrid(); - const int rLocalOffset = subGrid.getLocalDomain().offset[axis_element.space]; - const int rLocalSize = int(hBuffer.size().y() - 2 * rGuardCells); - const int rGlobalSize = subGrid.getGlobalDomain().size[axis_element.space]; - PMACC_VERIFY(rLocalSize == subGrid.getLocalDomain().size[axis_element.space]); - - /* globalDomain of the phase space */ - splash::Dimensions globalPhaseSpace_size(hBuffer.size().x(), rGlobalSize, 1); - - /* global moving window meta information */ - splash::Dimensions globalPhaseSpace_offset(0, 0, 0); - int globalMovingWindowOffset = 0; - int globalMovingWindowSize = rGlobalSize; - if(axis_element.space == AxisDescription::y) /* spatial axis == y */ - { - globalPhaseSpace_offset.set(0, numSlides * rLocalSize, 0); - Window window = MovingWindow::getInstance().getWindow(currentStep); - globalMovingWindowOffset = window.globalDimensions.offset[axis_element.space]; - globalMovingWindowSize = window.globalDimensions.size[axis_element.space]; - } - - /* localDomain: offset of it in the globalDomain and size */ - splash::Dimensions localPhaseSpace_offset(0, rLocalOffset, 0); - splash::Dimensions localPhaseSpace_size(hBuffer.size().x(), rLocalSize, 1); - - /** Dataset Name **************************************************/ - std::ostringstream dataSetName; - /* xpx or ypz or ... */ - dataSetName << fCoords.at(axis_element.space) << "p" << fCoords.at(axis_element.momentum); - - /** debug log *****************************************************/ - int rank; - MPI_CHECK(MPI_Comm_rank(mpiComm, &rank)); - log( - "Dump buffer %1% to %2% at offset %3% with size %4% for total size %5% for rank %6% / %7%") - % (*(hBuffer.origin()(0, rGuardCells))) % dataSetName.str() % localPhaseSpace_offset.toString() - % localPhaseSpace_size.toString() % globalPhaseSpace_size.toString() % rank % size; - - /** write local domain ********************************************/ - typename PICToSplash::type ctPhaseSpace; - - // avoid deadlock between not finished pmacc tasks and mpi calls in HDF5 - __getTransactionEvent().waitForFinished(); - - pdc.writeDomain( - currentStep, - /* global domain and my local offset within it */ - globalPhaseSpace_size, - localPhaseSpace_offset, - /* */ - ctPhaseSpace, - bufDim, - /* local data set dimensions */ - splash::Selection(localPhaseSpace_size), - /* data set name */ - dataSetName.str().c_str(), - /* global domain */ - splash::Domain(globalPhaseSpace_offset, globalPhaseSpace_size), - /* dataClass, buffer */ - DomainCollector::GridType, - &(*hBuffer.origin()(0, rGuardCells))); - - /** meta attributes for the data set: unit, range, moving window **/ - typedef PICToSplash::type SplashFloatXType; - typedef PICToSplash::type SplashFloat64Type; - ColTypeInt ctInt; - SplashFloat64Type ctFloat64; - SplashFloatXType ctFloatX; - - pdc.writeAttribute(currentStep, ctFloat64, dataSetName.str().c_str(), "sim_unit", &unit); - pdc.writeAttribute(currentStep, ctFloat64, dataSetName.str().c_str(), "p_unit", &pRange_unit); - pdc.writeAttribute(currentStep, ctFloatX, dataSetName.str().c_str(), "p_min", &(axis_p_range.first)); - pdc.writeAttribute(currentStep, ctFloatX, dataSetName.str().c_str(), "p_max", &(axis_p_range.second)); - pdc.writeAttribute( - currentStep, - ctInt, - dataSetName.str().c_str(), - "movingWindowOffset", - &globalMovingWindowOffset); - pdc.writeAttribute( - currentStep, - ctInt, - dataSetName.str().c_str(), - "movingWindowSize", - &globalMovingWindowSize); - - pdc.writeAttribute( - currentStep, - ctFloatX, - dataSetName.str().c_str(), - "dr", - &(cellSize[axis_element.space])); - pdc.writeAttribute(currentStep, ctFloatX, dataSetName.str().c_str(), "dV", &CELL_VOLUME); - pdc.writeAttribute(currentStep, ctFloat64, dataSetName.str().c_str(), "dr_unit", &UNIT_LENGTH); - pdc.writeAttribute(currentStep, ctFloatX, dataSetName.str().c_str(), "dt", &DELTA_T); - pdc.writeAttribute(currentStep, ctFloat64, dataSetName.str().c_str(), "dt_unit", &UNIT_TIME); - - /** close file ****************************************************/ - pdc.finalize(); - pdc.close(); - } - }; - -} /* namespace picongpu */ diff --git a/include/picongpu/plugins/PhaseSpace/PhaseSpace.hpp b/include/picongpu/plugins/PhaseSpace/PhaseSpace.hpp index e21eaf03108..438e71a5a3c 100644 --- a/include/picongpu/plugins/PhaseSpace/PhaseSpace.hpp +++ b/include/picongpu/plugins/PhaseSpace/PhaseSpace.hpp @@ -87,6 +87,19 @@ namespace picongpu plugins::multi::Option momentum_range_min = {"min", "min range momentum [m_species c]"}; plugins::multi::Option momentum_range_max = {"max", "max range momentum [m_species c]"}; + /* + * Set to h5 for now at least, to make for easier comparison of + * output with old outpu + */ + plugins::multi::Option file_name_extension + = {"ext", + "openPMD filename extension (this controls the" + "backend picked by the openPMD API)", + "h5"}; + + plugins::multi::Option json_config + = {"json", "advanced (backend) configuration for openPMD in JSON format", "{}"}; + //! string list with all possible particle filters std::string concatenatedFilterNames; std::vector allowedFilters; @@ -108,6 +121,8 @@ namespace picongpu element_momentum.registerHelp(desc, masterPrefix + prefix); momentum_range_min.registerHelp(desc, masterPrefix + prefix); momentum_range_max.registerHelp(desc, masterPrefix + prefix); + file_name_extension.registerHelp(desc, masterPrefix + prefix); + json_config.registerHelp(desc, masterPrefix + prefix); } void expandHelp( diff --git a/include/picongpu/plugins/PhaseSpace/PhaseSpace.tpp b/include/picongpu/plugins/PhaseSpace/PhaseSpace.tpp index 696fa1cf9d5..a92f5779d98 100644 --- a/include/picongpu/plugins/PhaseSpace/PhaseSpace.tpp +++ b/include/picongpu/plugins/PhaseSpace/PhaseSpace.tpp @@ -20,7 +20,7 @@ #pragma once #include "PhaseSpace.hpp" -#include "DumpHBufferSplashP.hpp" +#include "DumpHBufferOpenPMD.hpp" #include #include @@ -303,6 +303,8 @@ namespace picongpu pRange_unit, unit, Species::FrameType::getName() + "_" + m_help->filter.get(m_id), + m_help->file_name_extension.get(m_id), + m_help->json_config.get(m_id), currentStep, this->commFileWriter); } diff --git a/include/picongpu/plugins/PluginController.hpp b/include/picongpu/plugins/PluginController.hpp index 49ce302d258..42ef57928fa 100644 --- a/include/picongpu/plugins/PluginController.hpp +++ b/include/picongpu/plugins/PluginController.hpp @@ -211,6 +211,8 @@ namespace picongpu , plugins::radiation::Radiation, plugins::multi::Master>, +#endif +#if(ENABLE_OPENPMD == 1) plugins::multi::Master> #endif #if(PMACC_CUDA_ENABLED == 1) diff --git a/lib/python/picongpu/plugins/data/phase_space.py b/lib/python/picongpu/plugins/data/phase_space.py index 4109a35ccb0..5bad6fcf15e 100644 --- a/lib/python/picongpu/plugins/data/phase_space.py +++ b/lib/python/picongpu/plugins/data/phase_space.py @@ -10,9 +10,7 @@ import collections import numpy as np import os -import glob -import re -import h5py as h5 +import openpmd_api as io class PhaseSpaceMeta(object): @@ -71,10 +69,9 @@ def __init__(self, run_directory): super().__init__(run_directory) self.data_file_prefix = "PhaseSpace_{0}_{1}_{2}_{3}" - self.data_file_suffix = ".h5" self.data_hdf5_path = "/data/{0}/{1}" - def get_data_path(self, ps, species, species_filter="all", iteration=None): + def get_data_path(self, ps, species, species_filter="all", file_ext="h5"): """ Return the path to the underlying data file. @@ -89,19 +86,16 @@ def get_data_path(self, ps, species, species_filter="all", iteration=None): species_filter: string name of the particle species filter, default is 'all' (defined in ``particleFilters.param``) - iteration : (unsigned) int or list of int [unitless] - The iteration at which to read the data. - If 'None', a regular expression string matching - all iterations will be returned. + file_ext: string + filename extension for openPMD backend + default is 'h5' for the HDF5 backend Returns ------- - A string with a file path and a string with a in-file HDF5 path if - iteration is a single value or a list of length one. - If iteration is a list of length > 1, a list of paths is returned. - If iteration is None, only the first string is returned and contains a - regex-* for the position iteration. + A string with a the full openPMD file path pattern for loading from + a file-based iteration layout. """ + # @todo different file extensions? if species is None: raise ValueError('The species parameter can not be None!') if species_filter is None: @@ -123,45 +117,17 @@ def get_data_path(self, ps, species, species_filter="all", iteration=None): 'Did the simulation already run?' .format(self.run_directory)) - if iteration is not None: - if not isinstance(iteration, collections.Iterable): - iteration = [iteration] - - ret = [] - for it in iteration: - data_file_name = self.data_file_prefix.format( - species, - species_filter, - ps, - str(it)) + self.data_file_suffix - data_file_path = os.path.join(output_dir, data_file_name) - - if not os.path.isfile(data_file_path): - raise IOError('The file {} does not exist.\n' - 'Did the simulation already run?' - .format(data_file_path)) - - data_hdf5_name = self.data_hdf5_path.format( - it, - ps) - - ret.append((data_file_path, data_hdf5_name)) - if len(iteration) == 1: - return ret[0] - else: - return ret - else: - iteration_str = "*" - - data_file_name = self.data_file_prefix.format( - species, - species_filter, - ps, - iteration_str - ) + self.data_file_suffix - return os.path.join(output_dir, data_file_name) - - def get_iterations(self, ps, species, species_filter='all'): + iteration_str = "%T" + data_file_name = self.data_file_prefix.format( + species, + species_filter, + ps, + iteration_str + ) + '.' + file_ext + return os.path.join(output_dir, data_file_name) + + def get_iterations(self, ps, species, species_filter='all', + file_ext="h5"): """ Return an array of iterations with available data. @@ -176,32 +142,25 @@ def get_iterations(self, ps, species, species_filter='all'): species_filter: string name of the particle species filter, default is 'all' (defined in ``particleFilters.param``) + file_ext: string + filename extension for openPMD backend + default is 'h5' for the HDF5 backend Returns ------- An array with unsigned integers. """ # get the regular expression matching all available files - data_file_path = self.get_data_path(ps, species, species_filter) - - matching_files = glob.glob(data_file_path) - re_it = re.compile(data_file_path.replace("*", "([0-9]+)")) - - iterations = np.array( - sorted( - map( - lambda file_path: - np.uint64(re_it.match(file_path).group(1)), - matching_files - ) - ), - dtype=np.uint64 - ) + data_file_path = self.get_data_path(ps, species, species_filter, + file_ext=file_ext) + + series = io.Series(data_file_path, io.Access.read_only) + iterations = [key for key, _ in series.iterations.items()] return iterations def _get_for_iteration(self, iteration, ps, species, species_filter='all', - **kwargs): + file_ext="h5", **kwargs): """ Get a phase space histogram. @@ -219,6 +178,9 @@ def _get_for_iteration(self, iteration, ps, species, species_filter='all', species_filter: string name of the particle species filter, default is 'all' (defined in ``particleFilters.param``) + file_ext: string + filename extension for openPMD backend + default is 'h5' for the HDF5 backend Returns ------- @@ -231,8 +193,11 @@ def _get_for_iteration(self, iteration, ps, species, species_filter='all', containing ps and ps_meta for each requested iteration. If a single iteration is requested, return the tuple (ps, ps_meta). """ - available_iterations = self.get_iterations( - ps, species, species_filter) + + data_file_path = self.get_data_path(ps, species, species_filter, + file_ext=file_ext) + series = io.Series(data_file_path, io.Access.read_only) + available_iterations = [key for key, _ in series.iterations.items()] if iteration is not None: if not isinstance(iteration, collections.Iterable): @@ -247,28 +212,25 @@ def _get_for_iteration(self, iteration, ps, species, species_filter='all', iteration = available_iterations ret = [] - for it in iteration: - data_file_path, data_hdf5_name = self.get_data_path( - ps, - species, - species_filter, - it) - - f = h5.File(data_file_path, 'r') - ps_data = f[data_hdf5_name] + for index in iteration: + it = series.iterations[index] + dataset_name = "{}_{}_{}".format(species, species_filter, ps) + mesh = it.meshes[dataset_name] + ps_data = mesh[io.Mesh_Record_Component.SCALAR] # all in SI - dV = ps_data.attrs['dV'] * ps_data.attrs['dr_unit']**3 - unitSI = ps_data.attrs['sim_unit'] - p_range = ps_data.attrs['p_unit'] * \ - np.array([ps_data.attrs['p_min'], ps_data.attrs['p_max']]) - - mv_start = ps_data.attrs['movingWindowOffset'] - mv_end = mv_start + ps_data.attrs['movingWindowSize'] + dV = mesh.get_attribute('dV') * mesh.get_attribute('dr')**3 + unitSI = mesh.get_attribute('sim_unit') + p_range = mesh.get_attribute('p_unit') * \ + np.array( + [mesh.get_attribute('p_min'), mesh.get_attribute('p_max')]) + + mv_start = mesh.get_attribute('movingWindowOffset') + mv_end = mv_start + mesh.get_attribute('movingWindowSize') # 2D histogram: 0 (r_i); 1 (p_i) - spatial_offset = ps_data.attrs['_global_start'][1] + spatial_offset = mesh.get_attribute('_global_start')[0] - dr = ps_data.attrs['dr'] * ps_data.attrs['dr_unit'] + dr = mesh.get_attribute('dr') * mesh.get_attribute('dr_unit') r_range_cells = np.array([mv_start, mv_end]) + spatial_offset r_range = r_range_cells * dr @@ -278,7 +240,7 @@ def _get_for_iteration(self, iteration, ps, species, species_filter='all', # cut out the current window & scale by unitSI ps_cut = ps_data[mv_start:mv_end, :] * unitSI - f.close() + it.close() ps_meta = PhaseSpaceMeta( species, species_filter, ps, ps_cut.shape, extent, dV) diff --git a/lib/python/picongpu/plugins/plot_mpl/phase_space_visualizer.py b/lib/python/picongpu/plugins/plot_mpl/phase_space_visualizer.py index 85e6157e8fd..4c5cc163c1a 100644 --- a/lib/python/picongpu/plugins/plot_mpl/phase_space_visualizer.py +++ b/lib/python/picongpu/plugins/plot_mpl/phase_space_visualizer.py @@ -191,6 +191,9 @@ def visualize(self, **kwargs): ps : string phase space selection in order: spatial, momentum component, e.g. 'ypy' or 'ypx' + file_ext: string + filename extension for openPMD backend + default is 'h5' for the HDF5 backend """ super().visualize(**kwargs)