Main implementation

Use magic number instead of API call (impl)

Author: franzpoeschel

include/openPMD/Dataset.hpp

@@ -22,7 +22,9 @@
 
 #include "openPMD/Datatype.hpp"
 
+#include <limits>
 #include <memory>
+#include <optional>
 #include <string>
 #include <type_traits>
 #include <vector>
@@ -37,6 +39,11 @@ class Dataset
     friend class RecordComponent;
 
 public:
+    enum : std::uint64_t
+    {
+        JOINED_DIMENSION = std::numeric_limits<std::uint64_t>::max()
+    };
+
     Dataset(Datatype, Extent, std::string options = "{}");
 
     /**
@@ -53,5 +60,9 @@ class Dataset
     Datatype dtype;
     uint8_t rank;
     std::string options = "{}"; //!< backend-dependent JSON configuration
+
+    bool empty() const;
+
+    std::optional<size_t> joinedDimension() const;
 };
 } // namespace openPMD

include/openPMD/IO/IOTask.hpp

@@ -333,6 +333,7 @@ struct OPENPMDAPI_EXPORT Parameter<Operation::CREATE_DATASET>
     Extent extent = {};
     Datatype dtype = Datatype::UNDEFINED;
     std::string options = "{}";
+    std::optional<size_t> joinedDimension;
 
     /** Warn about unused JSON paramters
      *

include/openPMD/RecordComponent.tpp

@@ -259,8 +259,17 @@ RecordComponent::storeChunk(T_ContiguousContainer &data, Offset o, Extent e)
     // default arguments
     //   offset = {0u}: expand to right dim {0u, 0u, ...}
     Offset offset = o;
-    if (o.size() == 1u && o.at(0) == 0u && dim > 1u)
-        offset = Offset(dim, 0u);
+    if (o.size() == 1u && o.at(0) == 0u)
+    {
+        if (joinedDimension().has_value())
+        {
+            offset.clear();
+        }
+        else if (dim > 1u)
+        {
+            offset = Offset(dim, 0u);
+        }
+    }
 
     //   extent = {-1u}: take full size
     Extent extent(dim, 1u);
@@ -303,6 +312,7 @@ RecordComponent::storeChunk(Offset o, Extent e, F &&createBuffer)
         dCreate.name = rc.m_name;
         dCreate.extent = getExtent();
         dCreate.dtype = getDatatype();
+        dCreate.joinedDimension = joinedDimension();
         if (!rc.m_dataset.has_value())
         {
             throw error::WrongAPIUsage(

include/openPMD/backend/BaseRecordComponent.hpp

@@ -143,6 +143,8 @@ class BaseRecordComponent : virtual public Attributable
      */
     bool constant() const;
 
+    std::optional<size_t> joinedDimension() const;
+
     /**
      * Get data chunks that are available to be loaded from the backend.
      * Note that this is backend-dependent information and the returned

include/openPMD/backend/PatchRecordComponent.hpp

@@ -20,6 +20,7 @@
  */
 #pragma once
 
+#include "openPMD/Error.hpp"
 #include "openPMD/auxiliary/ShareRawInternal.hpp"
 #include "openPMD/backend/BaseRecordComponent.hpp"
 
@@ -115,6 +116,9 @@ class PatchRecordComponent : public BaseRecordComponent
     template <typename T>
     void store(uint64_t idx, T);
 
+    template <typename T>
+    void store(T);
+
     // clang-format off
 OPENPMD_private
     // clang-format on
@@ -244,4 +248,33 @@ inline void PatchRecordComponent::store(uint64_t idx, T data)
     auto &rc = get();
     rc.m_chunks.push(IOTask(this, std::move(dWrite)));
 }
+
+template <typename T>
+inline void PatchRecordComponent::store(T data)
+{
+    Datatype dtype = determineDatatype<T>();
+    if (dtype != getDatatype())
+    {
+        std::ostringstream oss;
+        oss << "Datatypes of patch data (" << dtype << ") and dataset ("
+            << getDatatype() << ") do not match.";
+        throw std::runtime_error(oss.str());
+    }
+
+    if (!joinedDimension().has_value())
+    {
+        throw error::WrongAPIUsage(
+            "[PatchRecordComponent::store] API call without explicit "
+            "specification of index only allowed when a joined dimension is "
+            "specified.");
+    }
+
+    Parameter<Operation::WRITE_DATASET> dWrite;
+    dWrite.offset = {};
+    dWrite.extent = {1};
+    dWrite.dtype = dtype;
+    dWrite.data = std::make_shared<T>(data);
+    auto &rc = get();
+    rc.m_chunks.push(IOTask(this, std::move(dWrite)));
+}
 } // namespace openPMD

src/Dataset.cpp

@@ -19,6 +19,7 @@
  * If not, see <http://www.gnu.org/licenses/>.
  */
 #include "openPMD/Dataset.hpp"
+#include "openPMD/Error.hpp"
 
 #include <cstddef>
 #include <iostream>
@@ -30,6 +31,9 @@ Dataset::Dataset(Datatype d, Extent e, std::string options_in)
 {
     // avoid initialization order issues
     rank = static_cast<uint8_t>(extent.size());
+    // Call this in order to have early error message in case of wrong
+    // specification of joined dimensions
+    joinedDimension();
 }
 
 Dataset::Dataset(Extent e) : Dataset(Datatype::UNDEFINED, std::move(e))
@@ -49,4 +53,41 @@ Dataset &Dataset::extend(Extent newExtents)
     extent = newExtents;
     return *this;
 }
+
+bool Dataset::empty() const
+{
+    auto jd = joinedDimension();
+    for (size_t i = 0; i < extent.size(); ++i)
+    {
+        if (extent[i] == 0 && (!jd.has_value() || jd.value() != i))
+        {
+            return true;
+        }
+    }
+    return false;
+}
+
+std::optional<size_t> Dataset::joinedDimension() const
+{
+    std::optional<size_t> res;
+    for (size_t i = 0; i < extent.size(); ++i)
+    {
+        if (extent[i] == JOINED_DIMENSION)
+        {
+            if (res.has_value())
+            {
+                throw error::WrongAPIUsage(
+                    "Must specify JOINED_DIMENSION at most once (found at "
+                    "indices " +
+                    std::to_string(res.value()) + " and " + std::to_string(i) +
+                    ")");
+            }
+            else
+            {
+                res = i;
+            }
+        }
+    }
+    return res;
+}
 } // namespace openPMD

src/IO/ADIOS/ADIOS2IOHandler.cpp

@@ -65,6 +65,21 @@ namespace openPMD
 
 #if openPMD_HAVE_ADIOS2
 
+namespace
+{
+    std::optional<size_t> joinedDimension(adios2::Dims const &dims)
+    {
+        for (size_t i = 0; i < dims.size(); ++i)
+        {
+            if (dims[i] == adios2::JoinedDim)
+            {
+                return i;
+            }
+        }
+        return std::nullopt;
+    }
+} // namespace
+
 #if openPMD_HAVE_MPI
 
 ADIOS2IOHandlerImpl::ADIOS2IOHandlerImpl(
@@ -720,8 +735,11 @@ void ADIOS2IOHandlerImpl::createDataset(
                 varName + "' remain unused:\n");
 
         // cast from openPMD::Extent to adios2::Dims
-        adios2::Dims const shape(
-            parameters.extent.begin(), parameters.extent.end());
+        adios2::Dims shape(parameters.extent.begin(), parameters.extent.end());
+        if (auto jd = parameters.joinedDimension; jd.has_value())
+        {
+            shape[jd.value()] = adios2::JoinedDim;
+        }
 
         auto &fileData = getFileData(file, IfFileNotOpen::ThrowError);
 
@@ -1623,11 +1641,22 @@ adios2::Variable<T> ADIOS2IOHandlerImpl::verifyDataset(
                 std::to_string(requiredDim) + ", but has dimensionality " +
                 std::to_string(actualDim) + ")")
     }
+    auto joinedDim = joinedDimension(shape);
     for (unsigned int i = 0; i < actualDim; i++)
+    {
+        if (!joinedDim.has_value() || i != joinedDim.value())
+        {
+            VERIFY_ALWAYS(
+                offset[i] + extent[i] <= shape[i],
+                "[ADIOS2] Dataset access out of bounds.")
+        }
+    }
+
+    if (joinedDim.has_value())
     {
         VERIFY_ALWAYS(
-            offset[i] + extent[i] <= shape[i],
-            "[ADIOS2] Dataset access out of bounds.")
+            offset.empty(),
+            "[ADIOS2] Offset must be an empty vector in case of joined array.");
     }
 
     var.SetSelection(

src/RecordComponent.cpp

@@ -95,10 +95,7 @@ RecordComponent &RecordComponent::resetDataset(Dataset d)
     }
     // if( d.extent.empty() )
     //    throw std::runtime_error("Dataset extent must be at least 1D.");
-    if (std::any_of(
-            d.extent.begin(), d.extent.end(), [](Extent::value_type const &i) {
-                return i == 0u;
-            }))
+    if (d.empty())
         return makeEmpty(std::move(d));
 
     rc.m_isEmpty = false;
@@ -299,6 +296,7 @@ void RecordComponent::flush(
                 dCreate.extent = getExtent();
                 dCreate.dtype = getDatatype();
                 dCreate.options = rc.m_dataset.value().options;
+                dCreate.joinedDimension = joinedDimension();
                 IOHandler()->enqueue(IOTask(this, dCreate));
             }
         }
@@ -474,22 +472,57 @@ void RecordComponent::verifyChunk(
     uint8_t dim = getDimensionality();
     Extent dse = getExtent();
 
-    if (e.size() != dim || o.size() != dim)
+    if (auto jd = joinedDimension(); jd.has_value())
     {
-        std::ostringstream oss;
-        oss << "Dimensionality of chunk ("
-            << "offset=" << o.size() << "D, "
-            << "extent=" << e.size() << "D) "
-            << "and record component (" << int(dim) << "D) "
-            << "do not match.";
-        throw std::runtime_error(oss.str());
+        if (o.size() != 0)
+        {
+            std::ostringstream oss;
+            oss << "Joined array: Must specify an empty offset (given: "
+                << "offset=" << o.size() << "D, "
+                << "extent=" << e.size() << "D).";
+            throw std::runtime_error(oss.str());
+        }
+        if (e.size() != dim)
+        {
+            std::ostringstream oss;
+            oss << "Joined array: Dimensionalities of chunk extent and dataset "
+                   "extent must be equivalent (given: "
+                << "offset=" << o.size() << "D, "
+                << "extent=" << e.size() << "D).";
+            throw std::runtime_error(oss.str());
+        }
+        for (size_t i = 0; i < dim; ++i)
+        {
+            if (i != jd.value() && e[i] != dse[i])
+            {
+                throw std::runtime_error(
+                    "Joined array: Chunk extent on non-joined dimensions must "
+                    "be equivalent to dataset extents (Dimension on index " +
+                    std::to_string(i) + ". DS: " + std::to_string(dse[i]) +
+                    " - Chunk: " + std::to_string(o[i] + e[i]) + ")");
+            }
+        }
+    }
+    else
+    {
+        if (e.size() != dim || o.size() != dim)
+        {
+            std::ostringstream oss;
+            oss << "Dimensionality of chunk ("
+                << "offset=" << o.size() << "D, "
+                << "extent=" << e.size() << "D) "
+                << "and record component (" << int(dim) << "D) "
+                << "do not match.";
+            throw std::runtime_error(oss.str());
+        }
+        for (uint8_t i = 0; i < dim; ++i)
+            if (dse[i] < o[i] + e[i])
+                throw std::runtime_error(
+                    "Chunk does not reside inside dataset (Dimension on "
+                    "index " +
+                    std::to_string(i) + ". DS: " + std::to_string(dse[i]) +
+                    " - Chunk: " + std::to_string(o[i] + e[i]) + ")");
     }
-    for (uint8_t i = 0; i < dim; ++i)
-        if (dse[i] < o[i] + e[i])
-            throw std::runtime_error(
-                "Chunk does not reside inside dataset (Dimension on index " +
-                std::to_string(i) + ". DS: " + std::to_string(dse[i]) +
-                " - Chunk: " + std::to_string(o[i] + e[i]) + ")");
 }
 
 namespace

src/backend/BaseRecordComponent.cpp

@@ -65,6 +65,19 @@ bool BaseRecordComponent::constant() const
     return get().m_isConstant;
 }
 
+std::optional<size_t> BaseRecordComponent::joinedDimension() const
+{
+    auto &rc = get();
+    if (rc.m_dataset.has_value())
+    {
+        return rc.m_dataset.value().joinedDimension();
+    }
+    else
+    {
+        return false;
+    }
+}
+
 ChunkTable BaseRecordComponent::availableChunks()
 {
     auto &rc = get();

src/backend/PatchRecordComponent.cpp

@@ -45,14 +45,11 @@ PatchRecordComponent &PatchRecordComponent::resetDataset(Dataset d)
             "written.");
     if (d.extent.empty())
         throw std::runtime_error("Dataset extent must be at least 1D.");
-    if (std::any_of(
-            d.extent.begin(), d.extent.end(), [](Extent::value_type const &i) {
-                return i == 0u;
-            }))
+    if (d.empty())
         throw std::runtime_error(
             "Dataset extent must not be zero in any dimension.");
 
-    get().m_dataset = d;
+    get().m_dataset = std::move(d);
     dirty() = true;
     return *this;
 }
@@ -136,6 +133,7 @@ void PatchRecordComponent::flush(
             dCreate.extent = getExtent();
             dCreate.dtype = getDatatype();
             dCreate.options = rc.m_dataset.value().options;
+            dCreate.joinedDimension = joinedDimension();
             IOHandler()->enqueue(IOTask(this, dCreate));
         }