Non-spatial meshes (#1534)

* Implement per-dimension gridUnitSI * Add gridUnitDimension * Add an alias API + some bit of testing * Add upgrade notice for old API call We cannot upgrade users silently in this place * Use openPMD 2.0 standard setting * Two little fixes * Default value for gridUnitSI set upon flush * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fallback for undefined gridUnitSI * Add conversion helpers between map- and array representations * Windows fixes * Tests and bugfixes * CI fixes * Introduce a member that stores the openPMD version as enum * Use enum-type standard check for meshes * Move auxiliary namespace under unit_representation * Remove retrieveSeries_optional * Cleaner Python API, Python tests * Python documentation, C++ documentation fixes * Better function argument types * CI fixes * Extend Python test --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
openPMD · Nov 12, 2024 · d1015ee · d1015ee
1 parent 4a040bb
commit d1015ee
Show file tree

Hide file tree

Showing 21 changed files with 839 additions and 88 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -401,6 +401,7 @@ set(CORE_SOURCE
         src/Record.cpp
         src/RecordComponent.cpp
         src/Series.cpp
+        src/UnitDimension.cpp
         src/version.cpp
         src/WriteIterations.cpp
         src/auxiliary/Date.cpp

diff --git a/include/openPMD/IO/AbstractIOHandler.hpp b/include/openPMD/IO/AbstractIOHandler.hpp
@@ -25,6 +25,7 @@
 #include "openPMD/IO/IOTask.hpp"
 #include "openPMD/IterationEncoding.hpp"
 #include "openPMD/config.hpp"
+#include "openPMD/version.hpp"
 
 #if openPMD_HAVE_MPI
 #include <mpi.h>
@@ -75,6 +76,20 @@ enum class FlushLevel
     CreateOrOpenFiles
 };
 
+enum class OpenpmdStandard
+{
+    v_1_0_0,
+    v_1_0_1,
+    v_1_1_0,
+    v_2_0_0
+};
+
+namespace auxiliary
+{
+    auto parseStandard(std::string const &) -> OpenpmdStandard;
+    auto formatStandard(OpenpmdStandard) -> char const *;
+} // namespace auxiliary
+
 namespace internal
 {
     /**
@@ -189,38 +204,7 @@ class AbstractIOHandler
     friend class detail::ADIOS2File;
 
 private:
-    IterationEncoding m_encoding = IterationEncoding::groupBased;
-
-    void setIterationEncoding(IterationEncoding encoding)
-    {
-        /*
-         * In file-based iteration encoding, the APPEND mode is handled entirely
-         * by the frontend, the backend should just treat it as CREATE mode.
-         * Similar for READ_LINEAR which should be treated as READ_RANDOM_ACCESS
-         * in the backend.
-         */
-        if (encoding == IterationEncoding::fileBased)
-        {
-            switch (m_backendAccess)
-            {
-
-            case Access::READ_LINEAR:
-                // do we really want to have those as const members..?
-                *const_cast<Access *>(&m_backendAccess) =
-                    Access::READ_RANDOM_ACCESS;
-                break;
-            case Access::APPEND:
-                *const_cast<Access *>(&m_backendAccess) = Access::CREATE;
-                break;
-            case Access::READ_RANDOM_ACCESS:
-            case Access::READ_WRITE:
-            case Access::CREATE:
-                break;
-            }
-        }
-
-        m_encoding = encoding;
-    }
+    void setIterationEncoding(IterationEncoding encoding);
 
 public:
 #if openPMD_HAVE_MPI
@@ -284,8 +268,14 @@ class AbstractIOHandler
      */
     Access m_backendAccess;
     Access m_frontendAccess;
-    internal::SeriesStatus m_seriesStatus = internal::SeriesStatus::Default;
     std::queue<IOTask> m_work;
+
+    /**************************************************************************
+     * Since the AbstractIOHandler is linked to every object of the frontend, *
+     * it stores a number of members that are needed by methods traversing    *
+     * the object hierarchy. Those members are found below.                   *
+     **************************************************************************/
+
     /**
      * This is to avoid that the destructor tries flushing again if an error
      * happened. Otherwise, this would lead to confusing error messages.
@@ -294,6 +284,9 @@ class AbstractIOHandler
      * The destructor will only attempt flushing again if this is true.
      */
     bool m_lastFlushSuccessful = false;
+    internal::SeriesStatus m_seriesStatus = internal::SeriesStatus::Default;
+    IterationEncoding m_encoding = IterationEncoding::groupBased;
+    OpenpmdStandard m_standard = auxiliary::parseStandard(getStandardDefault());
 }; // AbstractIOHandler
 
 } // namespace openPMD
diff --git a/include/openPMD/Mesh.hpp b/include/openPMD/Mesh.hpp
@@ -20,11 +20,11 @@
  */
 #pragma once
 
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/Attributable.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 #include "openPMD/backend/MeshRecordComponent.hpp"
 
-#include <array>
 #include <ostream>
 #include <string>
 #include <type_traits>
@@ -155,7 +155,7 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     template <
         typename T,
-        typename = std::enable_if_t<std::is_floating_point<T>::value> >
+        typename = std::enable_if_t<std::is_floating_point<T>::value>>
     Mesh &setGridSpacing(std::vector<T> const &gridSpacing);
 
     /**
@@ -184,22 +184,119 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
      *
+     * Valid for openPMD version 1.*.
+     * In order to specify the gridUnitSI per dimension (openPMD 2.*),
+     * use the vector overload or `setGridUnitSIPerDimension()`.
+     *
      * @param   gridUnitSI  unit-conversion factor to multiply each value in
      * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
      * simulation units to SI units.
      * @return  Reference to modified mesh.
      */
     Mesh &setGridUnitSI(double gridUnitSI);
 
+    /** Alias for `setGridUnitSIPerDimension(std::vector<double>)`.
+     *
+     * Set the unit-conversion factors per axis to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * Valid for openPMD 2.*.
+     * The legacy behavior (openPMD 1.*, a scalar gridUnitSI) is implemented
+     * by `setGridUnitSI(double)`.
+     *
+     * @param   gridUnitSI  unit-conversion factor to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * @return  Reference to modified mesh.
+     */
+    Mesh &setGridUnitSI(std::vector<double> const &gridUnitSI);
+
+    /**
+     * @return  A vector of the gridUnitSI per grid axis in the order of
+     * the axisLabels. If the gridUnitSI is defined as a scalar
+     * (legacy openPMD), the dimensionality is determined and a vector of
+     * `dimensionality` times the scalar vector is returned.
+     */
+    std::vector<double> gridUnitSIPerDimension() const;
+
+    /* Set the unit-conversion factors per axis to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * Valid for openPMD 2.*.
+     * The legacy behavior (openPMD 1.*, a scalar gridUnitSI) is implemented
+     * by `setGridUnitSI(double)`.
+     *
+     * @param   gridUnitSI  unit-conversion factor to multiply each value in
+     * Mesh::gridSpacing and Mesh::gridGlobalOffset, in order to convert from
+     * simulation units to SI units.
+     *
+     * @return  Reference to modified mesh.
+     */
+    Mesh &setGridUnitSIPerDimension(std::vector<double> const &gridUnitSI);
+
     /** Set the powers of the 7 base measures characterizing the record's unit
      * in SI.
      *
      * @param   unitDimension   map containing pairs of (UnitDimension, double)
      * that represent the power of the particular base.
      * @return  Reference to modified mesh.
      */
+    Mesh &setUnitDimension(unit_representations::AsMap const &unitDimension);
+
+    /** Set the powers of the 7 base measures characterizing the record's unit
+     * in SI.
+     *
+     * @param   unitDimension   array containing seven doubles, each
+     * representing the power of the particular base in order.
+     * @return  Reference to modified mesh.
+     */
+    Mesh &setUnitDimension(unit_representations::AsArray const &unitDimension);
+
+    /**
+     * @brief Set the unitDimension for each axis of the current grid.
+     *
+     * @param gridUnitDimension A vector of the unitDimensions for each
+     * axis of the grid in the order of the axisLabels, in dict representation.
+
+     * Behavior note: This is an updating method, meaning that an SI unit that
+     * has been defined before and is in the next call not explicitly set
+     * in the `std::map<UnitDimension, double>` will keep its previous value.
+     *
+     * @return Reference to modified mesh.
+     */
     Mesh &
-    setUnitDimension(std::map<UnitDimension, double> const &unitDimension);
+    setGridUnitDimension(unit_representations::AsMaps const &gridUnitDimension);
+
+    /**
+     * @brief Set the unitDimension for each axis of the current grid.
+     *
+     * @param gridUnitDimension A vector of the unitDimensions for each
+     * axis of the grid in the order of the axisLabels, in array representation.
+
+     * Behavior note: This is an updating method, meaning that an SI unit that
+     * has been defined before and is in the next call not explicitly set
+     * in the `std::map<UnitDimension, double>` will keep its previous value.
+     *
+     * @return Reference to modified mesh.
+     */
+    Mesh &setGridUnitDimension(
+        unit_representations::AsArrays const &gridUnitDimension);
+
+    /**
+     * @brief Return the physical dimensions of the mesh axes.
+
+     * If the attribute is not defined, the axes are assumed to be spatial
+     * and the return value will be according to this assumption.
+     * If the attribute is defined, the dimensionality of the return value is
+     * not checked against the dimensionality of the mesh.
+     *
+     * @return A vector of arrays, each array representing the SI unit of one
+     * mesh axis.
+     */
+    unit_representations::AsArrays gridUnitDimension() const;
 
     /**
      * @tparam  T   Floating point type of user-selected precision (e.g. float,
@@ -222,7 +319,7 @@ class Mesh : public BaseRecord<MeshRecordComponent>
      */
     template <
         typename T,
-        typename = std::enable_if_t<std::is_floating_point<T>::value> >
+        typename = std::enable_if_t<std::is_floating_point<T>::value>>
     Mesh &setTimeOffset(T timeOffset);
 
 private:

diff --git a/include/openPMD/Record.hpp b/include/openPMD/Record.hpp
@@ -21,9 +21,9 @@
 #pragma once
 
 #include "openPMD/RecordComponent.hpp"
+#include "openPMD/UnitDimension.hpp"
 #include "openPMD/backend/BaseRecord.hpp"
 
-#include <map>
 #include <string>
 #include <type_traits>
 
@@ -40,7 +40,8 @@ class Record : public BaseRecord<RecordComponent>
     Record &operator=(Record const &) = default;
     ~Record() override = default;
 
-    Record &setUnitDimension(std::map<UnitDimension, double> const &);
+    Record &setUnitDimension(unit_representations::AsMap const &);
+    Record &setUnitDimension(unit_representations::AsArray const &);
 
     template <typename T>
     T timeOffset() const;

diff --git a/include/openPMD/UnitDimension.hpp b/include/openPMD/UnitDimension.hpp
@@ -20,10 +20,16 @@
  */
 #pragma once
 
+#include <array>
 #include <cstdint>
+#include <map>
+#include <vector>
 
 namespace openPMD
 {
+
+using UnitDimensionExponent = double;
+
 /** Physical dimension of a record
  *
  * Dimensional base quantities of the international system of quantities
@@ -38,4 +44,25 @@ enum class UnitDimension : uint8_t
     N, //!< amount of substance
     J //!< luminous intensity
 };
+
+namespace unit_representations
+{
+    using AsMap = std::map<UnitDimension, UnitDimensionExponent>;
+    using AsArray = std::array<UnitDimensionExponent, 7>;
+
+    using AsMaps = std::vector<AsMap>;
+    using AsArrays = std::vector<AsArray>;
+
+    auto asArray(AsMap const &) -> AsArray;
+    auto asMap(AsArray const &, bool skip_zeros = true) -> AsMap;
+
+    auto asArrays(AsMaps const &) -> AsArrays;
+    auto asMaps(AsArrays const &, bool skip_zeros = true) -> AsMaps;
+
+    namespace auxiliary
+    {
+        void fromMapOfUnitDimension(
+            double *cursor, std::map<UnitDimension, double> const &udim);
+    } // namespace auxiliary
+} // namespace unit_representations
 } // namespace openPMD
diff --git a/include/openPMD/backend/Attributable.hpp b/include/openPMD/backend/Attributable.hpp
@@ -328,6 +328,8 @@ class Attributable
      */
     void touch();
 
+    [[nodiscard]] OpenpmdStandard openPMDStandard() const;
+
     // clang-format off
 OPENPMD_protected
     // clang-format on

diff --git a/include/openPMD/backend/BaseRecord.hpp b/include/openPMD/backend/BaseRecord.hpp
@@ -474,7 +474,7 @@ class BaseRecord
      *
      * @return powers of the 7 base measures in the order specified above
      */
-    std::array<double, 7> unitDimension() const;
+    unit_representations::AsArray unitDimension() const;
 
     void setDatasetDefined(BaseRecordComponent::Data_t &data) override
     {
@@ -928,7 +928,7 @@ auto BaseRecord<T_elem>::emplace(Args &&...args) -> std::pair<iterator, bool>
 }
 
 template <typename T_elem>
-inline std::array<double, 7> BaseRecord<T_elem>::unitDimension() const
+inline unit_representations::AsArray BaseRecord<T_elem>::unitDimension() const
 {
     return this->getAttribute("unitDimension")
         .template get<std::array<double, 7>>();

diff --git a/include/openPMD/binding/python/UnitDimension.hpp b/include/openPMD/binding/python/UnitDimension.hpp
@@ -24,7 +24,7 @@ namespace openPMD
 {
 namespace python
 {
-    constexpr auto doc_unit_dimension = R"docstr(
+    constexpr auto doc_unit_dimension = &R"docstr(
 Return the physical dimension (quantity) of a record
 
 Annotating the physical dimension of a record allows us to read data
@@ -40,7 +40,29 @@ See https://en.wikipedia.org/wiki/International_System_of_Quantities#Base_quanti
 See https://github.com/openPMD/openPMD-standard/blob/1.1.0/STANDARD.md#required-for-each-record
 
 Returns the powers of the 7 base measures in the order specified above.
-)docstr";
+)docstr"[1];
+
+    constexpr auto doc_mesh_unit_dimension = &R"docstr(
+Return the physical dimension (quantity) of the record axes
+
+Annotating the physical dimension of the record axes allows us to read data
+sets with arbitrary names and understand their purpose simply by
+dimensional analysis. The dimensional base quantities in openPMD are
+in order: length (L), mass (M), time (T), electric current (I),
+thermodynamic temperature (theta), amount of substance (N),
+luminous intensity (J) after the international system of quantities
+(ISQ).
+This attribute may be left out, the axes will then be interpreted as spatial.
+
+See https://en.wikipedia.org/wiki/Dimensional_analysis
+See https://en.wikipedia.org/wiki/International_System_of_Quantities#Base_quantities
+See https://github.com/openPMD/openPMD-standard/blob/1.1.0/STANDARD.md#required-for-each-record
+
+Returns the powers of the 7 base measures in the order specified above, listed
+for each axis in the order of the axisLabels.
+This attribute has been introduced as part of openPMD 2.0.0 in:
+Ref.: https://github.com/openPMD/openPMD-standard/pull/193
+)docstr"[1];
 
 } // namespace python
 } // namespace openPMD