Skip to content

Commit

Permalink
IO in PhaseSpace plugin via openPMD API
Browse files Browse the repository at this point in the history
Basic dataset dumping via openPMD

Fix dataset dimensions

const

Write some more attributes

Build PhaseSpace only if openPMD is available

Fix integer signedness

Use using

Create a new file for each species, enforce file-based iteration layout

Update naming once more

Make filenames dumped by openPMD identical to the ones formerly dumped
by splash. Give splash some legacy filename.

Move Python scripts to openPMD: get_iterations()

Move Python scripts to openPMD: _get_for_iteration

Cleanup python to make CI happy

Fix indexing

Update attributes naming

Formatting

Write _global_start and _global_offset

Go back to using scalar record components

Write some otherwise-defaulted openPMD attributes

Update documentation

Allow reading from any openPMD backend in the python scripts

Python formatting

Add (global|total)Domain(Size|Offset) Attributes

Write globalDomainAxisLabels

Document custom attributes and fix some errors

Remove libsplash output from phasespace plugin
  • Loading branch information
franzpoeschel committed Jan 27, 2021
1 parent a5582a0 commit 83570da
Show file tree
Hide file tree
Showing 9 changed files with 351 additions and 306 deletions.
26 changes: 20 additions & 6 deletions docs/source/usage/plugins/phaseSpace.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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 <install-dependencies>` are compiled in.
The plugin is available as soon as the :ref:`openPMD API <install-dependencies>` is compiled in.

.cfg file
^^^^^^^^^
Expand All @@ -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 <N>`` calculate each N steps *none*
``--e_phaseSpace.filter`` Use filtered particles. Available filters are set up in *none*
Expand All @@ -34,6 +34,7 @@ Option Usage
``--e_phaseSpace.momentum <px/py/pz>`` momentum coordinate of the 2D phase space *none*
``--e_phaseSpace.min <ValL>`` minimum of the momentum range :math:`m_\mathrm{species} c`
``--e_phaseSpace.max <ValR>`` maximum of the momentum range :math:`m_\mathrm{species} c`
``--e_phaseSpace.ext <ext>`` filename extension for openPMD backend *none*
====================================== ======================================================== ============================

Memory Complexity
Expand All @@ -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 <https://github.com/ComputationalRadiationPhysics/picongpu/wiki/PIConGPU-domain-definitions>`_ for explanations on the meaning of global and total domain.
* ``sim_unit``: SI scaling for the charge density values. Alias for ``unitSI``.

Analysis Tools
^^^^^^^^^^^^^^
Expand Down Expand Up @@ -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 <pp-openPMD>`.
- While the openPMD standard `has already been updated <https://github.com/openPMD/openPMD-standard/pull/193>`_ 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
^^^^^^^^^^
Expand Down
32 changes: 32 additions & 0 deletions include/picongpu/plugins/PhaseSpace/AxisDescription.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@

#pragma once

#include <string>

namespace picongpu
{
/** 2D Phase Space Selection
Expand Down Expand Up @@ -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 */
224 changes: 224 additions & 0 deletions include/picongpu/plugins/PhaseSpace/DumpHBufferOpenPMD.hpp
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "picongpu/simulation_defines.hpp"

#include "picongpu/plugins/PhaseSpace/AxisDescription.hpp"
#include <pmacc/communication/manager_common.hpp>
#include <pmacc/mappings/simulation/GridController.hpp>
#include <pmacc/mappings/simulation/SubGrid.hpp>
#include <pmacc/dimensions/DataSpace.hpp>
#include <pmacc/cuSTL/container/HostBuffer.hpp>
#include <pmacc/math/vector/Int.hpp>
#include <pmacc/verify.hpp>

#include <string>
#include <fstream>
#include <sstream>
#include <utility>
#include <mpi.h>
#include <openPMD/openPMD.hpp>
#include <vector>

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<typename T_Type, int T_bufDim>
void operator()(
const pmacc::container::HostBuffer<T_Type, T_bufDim>& hBuffer,
const AxisDescription axis_element,
const std::pair<float_X, float_X> 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<simDim>& gc = pmacc::Environment<simDim>::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<simDim>& subGrid = Environment<simDim>::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<picLog::INPUT_OUTPUT>(
"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<Type>(), globalPhaseSpace_extent});
std::shared_ptr<Type> data(&(*hBuffer.origin()(0, rGuardCells)), [](auto const&) {});
dataset.storeChunk<Type>(data, localPhaseSpace_offset, localPhaseSpace_extent);

/** meta attributes for the data set: unit, range, moving window **/

pmacc::Selection<simDim> globalDomain = subGrid.getGlobalDomain();
pmacc::Selection<simDim> totalDomain = subGrid.getTotalDomain();
// convert things to std::vector<> for the openPMD API to enjoy
std::vector<int> globalDomainSize{&globalDomain.size[0], &globalDomain.size[0] + simDim};
std::vector<int> globalDomainOffset{&globalDomain.offset[0], &globalDomain.offset[0] + simDim};
std::vector<int> totalDomainSize{&totalDomain.size[0], &totalDomain.size[0] + simDim};
std::vector<int> totalDomainOffset{&totalDomain.offset[0], &totalDomain.offset[0] + simDim};
std::vector<std::string> 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<float_X>{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<float>{0.5, 0.5});

/** close file ****************************************************/
iteration.close();
}
};

} /* namespace picongpu */
Loading

0 comments on commit 83570da

Please sign in to comment.