diff --git a/CMakeLists.txt b/CMakeLists.txt
index c95457bc6c19b6467eedabcb97f9d332d5442d1d..6727df3c90e685dd8146f68cebfaff32019bcf89 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -217,13 +217,33 @@ endif()
 
 if (${MPI_FOUND})
   set(PUGS_CXX_FLAGS "${PUGS_CXX_FLAGS} ${MPI_CXX_COMPILER_FLAGS}")
-  include_directories(SYSTEM ${MPI_CXX_INCLUDE_DIRS})
+  include_directories(SYSTEM "${MPI_CXX_INCLUDE_DIRS}")
 elseif(PUGS_ENABLE_MPI STREQUAL "ON")
   message(FATAL_ERROR "Could not find MPI library while requested")
 endif()
 
 set(PUGS_HAS_MPI ${MPI_FOUND})
 
+#------------------------------------------------------
+# search for HDF5
+
+set(PUGS_ENABLE_HDF5 AUTO CACHE STRING
+  "Choose one of: AUTO ON OFF")
+
+if (PUGS_ENABLE_HDF5 MATCHES "^(AUTO|ON)$")
+  set(HDF5_PREFER_PARALLEL true)
+  find_package(HDF5)
+
+  if (${HDF5_FOUND})
+    set(PUGS_CXX_FLAGS "${PUGS_CXX_FLAGS} ${HDF5_C_DEFINITIONS}")
+    include_directories(SYSTEM "${HDF5_C_INCLUDE_DIRS}")
+  endif()
+
+  set(PUGS_HAS_HDF5 ${HDF5_FOUND})
+else()
+  unset(PUGS_HAS_HDF5)
+endif()
+
 #------------------------------------------------------
 # search for clang-format
 
@@ -593,6 +613,8 @@ target_link_libraries(
   ${KOKKOS_CXX_FLAGS}
   ${OPENMP_LINK_FLAGS}
   ${PUGS_STD_LINK_FLAGS}
+  ${HDF5_C_LIBRARIES}
+  ${HDF5_C_HL_LIBRARIES}
   stdc++fs
   )
 
@@ -671,6 +693,16 @@ else()
   endif()
 endif()
 
+if (HDF5_FOUND)
+  message(" HDF5: ${HDF5_VERSION} parallel: ${HDF5_IS_PARALLEL}")
+else()
+  if (PUGS_ENABLE_SLEPC MATCHES "^(AUTO|ON)$")
+    message(" SLEPc: not found!")
+  else()
+      message(" SLEPc: explicitly deactivated!")
+  endif()
+endif()
+
 message("----------- utilities ----------")
 
 if(CLANG_FORMAT)
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 62dbd50536ad2e725e67a25d9cddf481a62b40c6..bf5eefe7cdbc21c463d962512a5ca00c8f0f846e 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -35,4 +35,5 @@ add_library(
   MeshNodeBoundary.cpp
   MeshRandomizer.cpp
   MeshTransformer.cpp
+  ParallelChecker.cpp
   SynchronizerManager.cpp)
diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 5a8fdb1bfdd66b209714a69b8f48c0de2fba07bd..5c922089e160bfd93d591bd99c153df9b9048ef5 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -17,6 +17,10 @@
 #include <slepc.h>
 #endif   // PUGS_HAS_PETSC
 
+#ifdef PUGS_HAS_HDF5
+#include <hdf5.h>
+#endif   // PUGS_HAS_HDF5
+
 std::string
 BuildInfo::type()
 {
@@ -73,3 +77,13 @@ BuildInfo::slepcLibrary()
   return "none";
 #endif   // PUGS_HAS_SLEPC
 }
+
+std::string
+BuildInfo::hdf5Library()
+{
+#ifdef PUGS_HAS_HDF5
+  return stringify(H5_VERSION) + ((H5_HAVE_PARALLEL) ? " [parallel]" : " [sequential]");
+#else
+  return "none";
+#endif   // PUGS_HAS_HDF5
+}
diff --git a/src/utils/BuildInfo.hpp b/src/utils/BuildInfo.hpp
index bc83cf3f1426bcab2972ebb8466f790cadf8a9ef..89a73e5411716c27a69de65fa9976fb062f77dca 100644
--- a/src/utils/BuildInfo.hpp
+++ b/src/utils/BuildInfo.hpp
@@ -11,6 +11,7 @@ struct BuildInfo
   static std::string mpiLibrary();
   static std::string petscLibrary();
   static std::string slepcLibrary();
+  static std::string hdf5Library();
 };
 
 #endif   // BUILD_INFO_HPP
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 2289a91d5f74c566a18af0f9fbcc8bf85ed844bb..7c6c99fcb8edee03b3cbed4687275eef3f98a527 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -8,6 +8,7 @@ add_library(
   Demangle.cpp
   Exceptions.cpp
   FPEManager.cpp
+  HDF5.cpp
   Messenger.cpp
   Partitioner.cpp
   PETScWrapper.cpp
diff --git a/src/utils/HDF5.cpp b/src/utils/HDF5.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4249d3520e3ce10a31b3419016d2bfd47288c088
--- /dev/null
+++ b/src/utils/HDF5.cpp
@@ -0,0 +1,187 @@
+#include <utils/HDF5.hpp>
+
+#ifdef PUGS_HAS_HDF5
+
+#include <utils/Exceptions.hpp>
+
+#include <rang.hpp>
+#include <sstream>
+
+HDF5::FileId
+HDF5::create(const std::string& filename)
+{
+  FileId file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+
+  if (file_id == -1) {
+    std::ostringstream error_msg;
+    error_msg << "Cannot create file '" << rang::fgB::yellow << filename << rang::fg::reset << "'";
+    throw NormalError(error_msg.str());
+  }
+  return file_id;
+}
+
+HDF5::FileId
+HDF5::openFileRW(const std::string& filename)
+{
+  FileId file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
+
+  if (file_id == -1) {
+    std::ostringstream error_msg;
+    error_msg << "Cannot open file '" << rang::fgB::yellow << filename << rang::fg::reset << "' for writing";
+    throw NormalError(error_msg.str());
+  }
+  return file_id;
+}
+
+HDF5::FileId
+HDF5::openFileRO(const std::string& filename)
+{
+  FileId file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+
+  if (file_id == -1) {
+    std::ostringstream error_msg;
+    error_msg << "Cannot read file '" << rang::fgB::yellow << filename << rang::fg::reset << "'";
+    throw NormalError(error_msg.str());
+  }
+  return file_id;
+}
+
+HDF5::GroupId
+HDF5::createOrOpenGroup(const FileId& file_id, const std::string& groupname)
+{
+  GroupId group_id = [&] {
+    if (H5Lexists(file_id, groupname.c_str(), H5P_DEFAULT) > 0) {
+      return H5Gopen2(file_id, groupname.c_str(), H5P_DEFAULT);
+    } else {
+      return H5Gcreate2(file_id, groupname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    }
+  }();
+
+  return group_id;
+}
+
+HDF5::GroupId
+HDF5::createOrOpenSubGroup(const GroupId& group_id, const std::string& groupname)
+{
+  GroupId sub_group_id = [&] {
+    if (H5Lexists(group_id, groupname.c_str(), H5P_DEFAULT) > 0) {
+      return H5Gopen2(group_id, groupname.c_str(), H5P_DEFAULT);
+    } else {
+      return H5Gcreate2(group_id, groupname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    }
+  }();
+
+  return sub_group_id;
+}
+
+HDF5::GroupId
+HDF5::openGroup(const FileId& file_id, const std::string& groupname)
+{
+  if (H5Lexists(file_id, groupname.c_str(), H5P_DEFAULT) > 0) {
+    return H5Gopen2(file_id, groupname.c_str(), H5P_DEFAULT);
+  } else {
+    std::ostringstream error_msg;
+    error_msg << "Cannot read group '" << rang::fgB::yellow << groupname << rang::fg::reset << "'";
+    throw NormalError(error_msg.str());
+  }
+}
+
+HDF5::GroupId
+HDF5::openSubGroup(const GroupId& group_id, const std::string& groupname)
+{
+  if (H5Lexists(group_id, groupname.c_str(), H5P_DEFAULT) > 0) {
+    return H5Gopen2(group_id, groupname.c_str(), H5P_DEFAULT);
+  } else {
+    std::ostringstream error_msg;
+    error_msg << "Cannot read sub-group '" << rang::fgB::yellow << groupname << rang::fg::reset << "'";
+    throw NormalError(error_msg.str());
+  }
+}
+
+template <HDF5::IdType id_type>
+void
+HDF5::_writeDataset(const HDFId<id_type>& id,
+                    const std::string& name,
+                    const DatatypeId& type,
+                    const size_t& array_size,
+                    const void* array)
+{
+  hsize_t dataspace_size[] = {array_size};
+  hid_t dataspace          = H5Screate_simple(1, dataspace_size, NULL);
+
+  hid_t dataset = H5Dcreate2(id, name.c_str(), type, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+  if (dataset == -1) {
+    throw NormalError("Cannot create dataset");
+  }
+  herr_t status = H5Dwrite(dataset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, array);
+  if (status < 0) {
+    throw NormalError("Cannot write dataset");
+  }
+
+  H5Dclose(dataset);
+  H5Sclose(dataspace);
+}
+
+template void HDF5::_writeDataset(const FileId& id,
+                                  const std::string& name,
+                                  const DatatypeId& type,
+                                  const size_t& array_size,
+                                  const void* array);
+
+template void HDF5::_writeDataset(const GroupId& id,
+                                  const std::string& name,
+                                  const DatatypeId& type,
+                                  const size_t& array_size,
+                                  const void* array);
+
+template <HDF5::IdType id_type>
+void
+HDF5::_writeAttribute(const HDFId<id_type>& id, const std::string& name, const DatatypeId& type, const void* value)
+{
+  hid_t space = H5Screate(H5S_SCALAR);
+
+  hid_t attribute = H5Acreate2(id, name.c_str(), type, space, H5P_DEFAULT, H5P_DEFAULT);
+
+  H5Awrite(attribute, type, value);
+}
+
+template void HDF5::_writeAttribute(const FileId& id,
+                                    const std::string& name,
+                                    const DatatypeId& type,
+                                    const void* value);
+
+template void HDF5::_writeAttribute(const GroupId& id,
+                                    const std::string& name,
+                                    const DatatypeId& type,
+                                    const void* value);
+
+template <HDF5::IdType id_type>
+void
+HDF5::_readAttribute(const HDFId<id_type>& id, const std::string& name, const DatatypeId& type, void* value)
+{
+  if (H5Aexists(id, name.c_str())) {
+    hid_t attribute_id   = H5Aopen(id, name.c_str(), H5P_DEFAULT);
+    hid_t attribute_type = H5Aget_type(attribute_id);
+    if (not H5Tequal(attribute_type, type)) {
+      std::ostringstream error_msg;
+      error_msg << "invalid type for attribute '" << rang::fgB::yellow << name << rang::fg::reset << "' got "
+                << H5Aget_type(attribute_id) << " expecting " << type;
+      throw NormalError(error_msg.str());
+    }
+
+    H5Tclose(attribute_type);
+    H5Aread(attribute_id, type, value);
+    H5Aclose(attribute_id);
+  } else {
+    std::ostringstream error_msg;
+    error_msg << "cannot find attribute '" << rang::fgB::yellow << name << rang::fg::reset << "'";
+    throw NormalError(error_msg.str());
+  }
+}
+
+template void HDF5::_readAttribute(const FileId& id, const std::string& name, const DatatypeId& type, void* value);
+
+template void HDF5::_readAttribute(const GroupId& id, const std::string& name, const DatatypeId& type, void* value);
+
+#endif   // PUGS_HAS_HDF5
diff --git a/src/utils/HDF5.hpp b/src/utils/HDF5.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c0e03503aa42c743e78c5cf85023c4000df548e8
--- /dev/null
+++ b/src/utils/HDF5.hpp
@@ -0,0 +1,308 @@
+#ifndef HDF5_HPP
+#define HDF5_HPP
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_HDF5
+
+#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <hdf5.h>
+
+#include <utils/Demangle.hpp>
+
+class HDF5
+{
+ private:
+  enum class IdType
+  {
+    attribute,
+    dataset,
+    file,
+    group,
+    datatype
+  };
+
+  template <IdType type>
+  class HDFId
+  {
+   private:
+    hid_t m_hdf_id;
+
+   public:
+    operator hid_t&()
+    {
+      return m_hdf_id;
+    }
+
+    operator hid_t() const
+    {
+      return m_hdf_id;
+    }
+
+    HDFId(hid_t hdf_id) : m_hdf_id{hdf_id} {}
+
+    template <IdType other_type>
+    HDFId(const HDFId<other_type>&)
+    {
+      static_assert(other_type == type, "HDF types conversion is forbidden");
+    }
+
+    template <IdType other_type>
+    HDFId(HDFId<other_type>&&)
+    {
+      static_assert(other_type == type, "HDF types conversion is forbidden");
+    }
+
+    HDFId(const HDFId&) = default;
+    HDFId(HDFId&&)      = default;
+
+    ~HDFId() = default;
+  };
+
+ public:
+  using AttributeId = HDFId<IdType::attribute>;
+  using DatasetId   = HDFId<IdType::dataset>;
+  using FileId      = HDFId<IdType::file>;
+  using GroupId     = HDFId<IdType::group>;
+  using DatatypeId  = HDFId<IdType::datatype>;
+
+ private:
+  template <IdType id_type>
+  static void _writeDataset(const HDFId<id_type>& id,
+                            const std::string& data_name,
+                            const DatatypeId& data_type,
+                            const size_t& array_size,
+                            const void* array);
+
+  template <IdType id_type>
+  static void _writeAttribute(const HDFId<id_type>& id,
+                              const std::string& data_name,
+                              const DatatypeId& data_type,
+                              const void* value);
+
+  template <IdType id_type>
+  static void _readAttribute(const HDFId<id_type>& id,
+                             const std::string& data_name,
+                             const DatatypeId& data_type,
+                             void* value);
+
+  template <typename DataT>
+  static DatatypeId
+  _getDataType()
+  {
+    using T = std::decay_t<DataT>;
+
+    static_assert(not std::is_pointer_v<DataT>);
+
+    if constexpr (std::is_same_v<char, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_CHAR);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<unsigned char, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_UCHAR);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<short, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_SHORT);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<unsigned short, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_USHORT);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<int, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_INT);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<unsigned int, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_UINT);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<long int, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_LONG);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<long unsigned int, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_ULONG);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<int8_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_INT8);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<uint8_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_UINT8);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<int16_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_INT16);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<uint16_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_UINT16);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<int32_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_INT32);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<uint32_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_UINT32);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<int64_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_INT64);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<uint64_t, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_UINT64);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+      //
+    } else if constexpr (std::is_same_v<float, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_FLOAT);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<double, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else if constexpr (std::is_same_v<long double, T>) {
+      DatatypeId datatype = H5Tcopy(H5T_NATIVE_LDOUBLE);
+      H5Tset_order(datatype, H5T_ORDER_LE);
+      return datatype;
+    } else {
+      //
+      static_assert(std::is_same_v<int, T>, "unexpected data type");
+    }
+  }
+
+ public:
+  static FileId create(const std::string& file_name);
+  static FileId openFileRW(const std::string& file_name);
+  static FileId openFileRO(const std::string& file_name);
+
+  static GroupId createOrOpenGroup(const FileId& file_id, const std::string& group_name);
+  static GroupId createOrOpenSubGroup(const GroupId& group_id, const std::string& group_name);
+
+  static GroupId openGroup(const FileId& file_id, const std::string& group_name);
+  static GroupId openSubGroup(const GroupId& group_id, const std::string& group_name);
+
+  template <typename DataT, IdType id_type>
+  static void
+  write(const HDFId<id_type>& id, const std::string& name, const Array<DataT>& array)
+  {
+    static_assert((id_type == IdType::file) or (id_type == IdType::group));
+    DatatypeId datatype = _getDataType<DataT>();
+    H5Tset_order(datatype, H5T_ORDER_LE);
+
+    const void* array_start = &(array[typename Array<DataT>::index_type(0)]);
+
+    _writeDataset(id, name, datatype, array.size(), array_start);
+
+    H5Tclose(datatype);
+  }
+
+  template <typename DataT, IdType id_type>
+  static Array<DataT>
+  readArray(const HDFId<id_type>& id, const std::string& name)
+  {
+    static_assert((id_type == IdType::file) or (id_type == IdType::group));
+    static_assert(not std::is_const_v<DataT>, "DataT must be mutable");
+
+    DatatypeId datatype = _getDataType<DataT>();
+    H5Tset_order(datatype, H5T_ORDER_LE);
+
+    DatasetId dataset = H5Dopen2(id, name.c_str(), H5P_DEFAULT);
+
+    if (dataset < 0) {
+      std::ostringstream error_msg;
+      error_msg << "cannot open dataset \"" << name << "\"\n";
+      throw NormalError(error_msg.str());
+    }
+
+    const size_t array_size = H5Dget_storage_size(dataset) / sizeof(DataT);
+    Array<DataT> array{array_size};
+
+    H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(array[0]));
+
+    H5Dclose(dataset);
+
+    H5Tclose(datatype);
+
+    return array;
+  }
+
+  template <typename DataT, IdType id_type>
+  static void
+  writeAttribute(const HDFId<id_type>& id, const std::string& name, const DataT& value)
+  {
+    static_assert((id_type == IdType::file) or (id_type == IdType::group));
+    if constexpr (std::is_same_v<std::string, std::decay_t<DataT>>) {
+      DatatypeId datatype = H5Tcopy(H5T_C_S1);
+      H5Tset_size(datatype, H5T_VARIABLE);
+
+      const void* v = value.c_str();
+      _writeAttribute(id, name, datatype, &v);
+      H5Tclose(datatype);
+    } else if constexpr (std::is_same_v<std::string_view, std::decay_t<DataT>>) {
+      DatatypeId datatype = H5Tcopy(H5T_C_S1);
+      H5Tset_size(datatype, H5T_VARIABLE);
+
+      const void* v = &(value[0]);
+      _writeAttribute(id, name, datatype, &v);
+      H5Tclose(datatype);
+    } else if constexpr (std::is_same_v<DataT, const char*> or std::is_same_v<DataT, char*>) {
+      DatatypeId datatype = H5Tcopy(H5T_C_S1);
+      H5Tset_size(datatype, H5T_VARIABLE);
+
+      _writeAttribute(id, name, datatype, &value);
+      H5Tclose(datatype);
+    } else if constexpr (std::is_arithmetic_v<std::decay_t<DataT>>) {
+      DatatypeId datatype = _getDataType<DataT>();
+      _writeAttribute(id, name, datatype, &value);
+      H5Tclose(datatype);
+    } else {
+      static_assert(std::is_same_v<std::decay_t<DataT>, std::string>, "unexpected attribute type");
+    }
+  }
+
+  template <typename DataT, IdType id_type>
+  static DataT
+  readAttribute(const HDFId<id_type>& id, const std::string& name)
+  {
+    static_assert(not std::is_const_v<DataT>, "cannot read to const value");
+    if constexpr (std::is_same_v<std::string, std::decay_t<DataT>>) {
+      DatatypeId datatype = H5Tcopy(H5T_C_S1);
+      H5Tset_size(datatype, H5T_VARIABLE);
+
+      char* v = nullptr;
+      _readAttribute(id, name, datatype, &v);
+
+      DataT value = v;
+      H5free_memory(v);
+
+      H5Tclose(datatype);
+      return value;
+    } else if constexpr (std::is_arithmetic_v<std::decay_t<DataT>>) {
+      DatatypeId datatype = _getDataType<DataT>();
+      DataT value;
+      _readAttribute(id, name, datatype, &value);
+      H5Tclose(datatype);
+      return value;
+    } else {
+      static_assert(std::is_same_v<std::decay_t<DataT>, std::string>, "unexpected attribute type");
+    }
+  }
+
+  static void
+  close(const FileId& file_id)
+  {
+    H5Fclose(file_id);
+  }
+};
+#endif   // PUGS_HAS_HDF5
+
+#endif   // HDF5_HPP
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index 4e1f521656683b371f488ee90df5d49e061d9c7b..5fcc6e8c48143496509421c3a0c8519e0dd23eea 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -59,6 +59,7 @@ pugsBuildInfo()
   os << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
   os << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
   os << "SLEPc:    " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
+  os << "HDF5:     " << rang::style::bold << BuildInfo::hdf5Library() << rang::style::reset << '\n';
   os << "-------------------------------------------------------";
 
   return os.str();
diff --git a/src/utils/pugs_config.hpp.in b/src/utils/pugs_config.hpp.in
index 0736003b9f40c4f71feea089072c19d9238e4033..4dc5babd4c157b9b41161dfbd3bb1156c1c6b631 100644
--- a/src/utils/pugs_config.hpp.in
+++ b/src/utils/pugs_config.hpp.in
@@ -5,6 +5,7 @@
 #cmakedefine PUGS_HAS_MPI
 #cmakedefine PUGS_HAS_PETSC
 #cmakedefine PUGS_HAS_SLEPC
+#cmakedefine PUGS_HAS_HDF5
 
 #cmakedefine SYSTEM_IS_LINUX
 #cmakedefine SYSTEM_IS_DARWIN
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9f1f7d9005d6aa8a583283bf2bfda8dc6e23eedd..0e2f09ffe07b34912dbef6e2715e7ee5e02f838e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -96,6 +96,7 @@ add_executable (unit_tests
   test_GaussLegendreQuadratureDescriptor.cpp
   test_GaussLobattoQuadratureDescriptor.cpp
   test_GaussQuadratureDescriptor.cpp
+  test_HDF5.cpp
   test_IfProcessor.cpp
   test_IncDecExpressionProcessor.cpp
   test_IntegrateCellArray.cpp
@@ -228,6 +229,8 @@ target_link_libraries (unit_tests
   ${PETSC_LIBRARIES}
   Catch2
   ${PUGS_STD_LINK_FLAGS}
+  ${HDF5_C_LIBRARIES}
+  ${HDF5_C_HL_LIBRARIES}
   stdc++fs
   )
 
@@ -255,6 +258,8 @@ target_link_libraries (mpi_unit_tests
   ${PETSC_LIBRARIES}
   Catch2
   ${PUGS_STD_LINK_FLAGS}
+  ${HDF5_C_LIBRARIES}
+  ${HDF5_C_HL_LIBRARIES}
   stdc++fs
   )
 
diff --git a/tests/test_HDF5.cpp b/tests/test_HDF5.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f84b5f694b84770ac8ba0406ed423f265f429f01
--- /dev/null
+++ b/tests/test_HDF5.cpp
@@ -0,0 +1,93 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/HDF5.hpp>
+
+#ifdef PUGS_HAS_HDF5
+
+#include <filesystem>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("HDF5", "[utils]")
+{
+  static const std::string filename = [&]() -> std::filesystem::path {
+    std::string template_temp_dir = std::filesystem::temp_directory_path() / "test_hdf5_XXXXXX";
+    return std::string{std::filesystem::path{mkdtemp(&template_temp_dir[0])}} + std::string{".h5"};
+  }();
+
+  SECTION("create file")
+  {
+    std::clog << "filename = " << filename << '\n';
+    HDF5::FileId file_id = HDF5::create(filename);
+    REQUIRE(file_id >= 0);
+    HDF5::close(file_id);
+  }
+
+  SECTION("open file and add attributes")
+  {
+    std::clog << "filename = " << filename << '\n';
+    HDF5::FileId file_id = HDF5::openFileRW(filename);
+    REQUIRE(file_id >= 0);
+
+    HDF5::writeAttribute(file_id, "char attr", static_cast<char>(-'b'));
+    HDF5::writeAttribute(file_id, "uchar attr", static_cast<unsigned char>('a'));
+    HDF5::writeAttribute(file_id, "short attr", static_cast<short>(-13));
+    HDF5::writeAttribute(file_id, "unsigned short attr", static_cast<unsigned short>(3));
+    HDF5::writeAttribute(file_id, "int attr", static_cast<int>(-11));
+    HDF5::writeAttribute(file_id, "unsigned int attr", static_cast<unsigned int>(23));
+    HDF5::writeAttribute(file_id, "long int attr", static_cast<long int>(-101));
+    HDF5::writeAttribute(file_id, "long unsigned int attr", static_cast<unsigned long int>(101));
+    HDF5::writeAttribute(file_id, "int8_t attr", static_cast<int8_t>(-25));
+    HDF5::writeAttribute(file_id, "uint8_t attr", static_cast<uint8_t>(27));
+    HDF5::writeAttribute(file_id, "int16_t attr", static_cast<int16_t>(-14));
+    HDF5::writeAttribute(file_id, "uint16_t attr", static_cast<uint16_t>(24));
+    HDF5::writeAttribute(file_id, "int32_t attr", static_cast<int32_t>(-124));
+    HDF5::writeAttribute(file_id, "uint32_t attr", static_cast<uint32_t>(147));
+    HDF5::writeAttribute(file_id, "int64_t attr", static_cast<int64_t>(-1515));
+    HDF5::writeAttribute(file_id, "uint64_t attr", static_cast<uint64_t>(1999));
+    HDF5::writeAttribute(file_id, "float attr", static_cast<float>(3.14));
+    HDF5::writeAttribute(file_id, "double attr", static_cast<double>(3.14159));
+    HDF5::writeAttribute(file_id, "long double attr", static_cast<long double>(3.14159265));
+    HDF5::writeAttribute(file_id, "foo attr", std::string{"foobar"});
+
+    HDF5::close(file_id);
+  }
+
+  SECTION("open file ro")
+  {
+    std::clog << "filename = " << filename << '\n';
+    HDF5::FileId file_id = HDF5::openFileRW(filename);
+    REQUIRE(file_id >= 0);
+
+    REQUIRE(HDF5::readAttribute<char>(file_id, "char attr") == static_cast<char>(-'b'));
+    REQUIRE(HDF5::readAttribute<unsigned char>(file_id, "uchar attr") == static_cast<unsigned char>('a'));
+    REQUIRE(HDF5::readAttribute<short>(file_id, "short attr") == static_cast<short>(-13));
+    REQUIRE(HDF5::readAttribute<unsigned short>(file_id, "unsigned short attr") == static_cast<unsigned short>(3));
+    REQUIRE(HDF5::readAttribute<int>(file_id, "int attr") == static_cast<int>(-11));
+    REQUIRE(HDF5::readAttribute<unsigned int>(file_id, "unsigned int attr") == static_cast<unsigned int>(23));
+    REQUIRE(HDF5::readAttribute<long int>(file_id, "long int attr") == static_cast<long int>(-101));
+    REQUIRE(HDF5::readAttribute<unsigned long int>(file_id, "long unsigned int attr") ==
+            static_cast<unsigned long int>(101));
+    REQUIRE(HDF5::readAttribute<int8_t>(file_id, "int8_t attr") == static_cast<int8_t>(-25));
+    REQUIRE(HDF5::readAttribute<uint8_t>(file_id, "uint8_t attr") == static_cast<uint8_t>(27));
+    REQUIRE(HDF5::readAttribute<int16_t>(file_id, "int16_t attr") == static_cast<int16_t>(-14));
+    REQUIRE(HDF5::readAttribute<uint16_t>(file_id, "uint16_t attr") == static_cast<uint16_t>(24));
+    REQUIRE(HDF5::readAttribute<int32_t>(file_id, "int32_t attr") == static_cast<int32_t>(-124));
+    REQUIRE(HDF5::readAttribute<uint32_t>(file_id, "uint32_t attr") == static_cast<uint32_t>(147));
+    REQUIRE(HDF5::readAttribute<int64_t>(file_id, "int64_t attr") == static_cast<int64_t>(-1515));
+    REQUIRE(HDF5::readAttribute<uint64_t>(file_id, "uint64_t attr") == static_cast<uint64_t>(1999));
+    REQUIRE(HDF5::readAttribute<float>(file_id, "float attr") == static_cast<float>(3.14));
+    REQUIRE(HDF5::readAttribute<double>(file_id, "double attr") == static_cast<double>(3.14159));
+    REQUIRE(HDF5::readAttribute<long double>(file_id, "long double attr") == static_cast<long double>(3.14159265));
+    REQUIRE(HDF5::readAttribute<std::string>(file_id, "foo attr") == std::string{"foobar"});
+
+    REQUIRE_THROWS_WITH(HDF5::readAttribute<double>(file_id, "invalid attr"),
+                        "error: cannot find attribute 'invalid attr'");
+    REQUIRE_THROWS_WITH(HDF5::readAttribute<double>(file_id, "int64_t attr"),
+                        "error: cannot find attribute 'invalid attr'");
+
+    HDF5::close(file_id);
+  }
+}
+#endif   // PUGS_HAS_HDF5
diff --git a/tests/test_PugsUtils.cpp b/tests/test_PugsUtils.cpp
index 19d693c76599937de4b789c2c07785844daa100c..9dfd686da65c5027f43bdfcf38ba3b54f90607ba 100644
--- a/tests/test_PugsUtils.cpp
+++ b/tests/test_PugsUtils.cpp
@@ -51,6 +51,7 @@ TEST_CASE("PugsUtils", "[utils]")
       os << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
       os << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
       os << "SLEPc:    " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
+      os << "HDF5:     " << rang::style::bold << BuildInfo::hdf5Library() << rang::style::reset << '\n';
       os << "-------------------------------------------------------";
 
       return os.str();