From b28ae4ad69a15bd063b8bd060c2582a1512f3508 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 19 Sep 2023 08:01:08 +0200
Subject: [PATCH] Switch to HighFive as an interface to HDF5

- add support for ItemArray
- began few code redesign/optimizations
---
 CMakeLists.txt                          |  63 +-
 src/CMakeLists.txt                      |   3 +
 src/dev_utils/CMakeLists.txt            |  11 +
 src/dev_utils/ParallelChecker.cpp       |  35 ++
 src/dev_utils/ParallelChecker.hpp       | 741 ++++++++++++++++++++++++
 src/language/modules/CMakeLists.txt     |  13 +-
 src/language/modules/DevUtilsModule.cpp |  35 ++
 src/language/utils/CMakeLists.txt       |  15 +-
 src/main.cpp                            |   6 +-
 src/mesh/CMakeLists.txt                 |   9 +-
 src/mesh/ParallelChecker.cpp            |  21 -
 src/mesh/ParallelChecker.hpp            | 363 ------------
 src/output/CMakeLists.txt               |   5 +
 src/scheme/CMakeLists.txt               |   8 +-
 src/scheme/DiscreteFunctionVariant.hpp  |   2 +-
 src/utils/BuildInfo.cpp                 |   2 +-
 src/utils/CMakeLists.txt                |   2 +-
 src/utils/HDF5.cpp                      | 235 --------
 src/utils/HDF5.hpp                      | 436 --------------
 tests/CMakeLists.txt                    |   9 +-
 tests/test_HDF5.cpp                     | 388 -------------
 21 files changed, 898 insertions(+), 1504 deletions(-)
 create mode 100644 src/dev_utils/CMakeLists.txt
 create mode 100644 src/dev_utils/ParallelChecker.cpp
 create mode 100644 src/dev_utils/ParallelChecker.hpp
 delete mode 100644 src/mesh/ParallelChecker.cpp
 delete mode 100644 src/mesh/ParallelChecker.hpp
 delete mode 100644 src/utils/HDF5.cpp
 delete mode 100644 src/utils/HDF5.hpp
 delete mode 100644 tests/test_HDF5.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1e0432ba3..ef0a0f4fa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -105,12 +105,6 @@ elseif(CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
   set(PUGS_CXX_FLAGS "${PUGS_CXX_FLAGS} -Wsign-compare -Wunused -Wunused-member-function -Wunused-private-field")
 endif()
 
-# if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
-#   if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "9.0.0")
-#     set(PUGS_STD_LINK_FLAGS "-lstdc++fs")
-#   endif()
-# endif()
-
 #------------------------------------------------------
 
 include (TestBigEndian)
@@ -178,7 +172,7 @@ else()
 endif()
 
 if (${PETSC_FOUND})
-  include_directories(SYSTEM ${PETSC_INCLUDE_DIRS})
+  include_directories(SYSTEM "${PETSC_INCLUDE_DIRS}")
 else()
   if (PUGS_ENABLE_PETSC MATCHES "^ON$")
     message(FATAL_ERROR "Could not find PETSc!")
@@ -206,7 +200,7 @@ else()
 endif()
 
 if (${SLEPC_FOUND})
-  include_directories(SYSTEM ${SLEPC_INCLUDE_DIRS})
+  include_directories(SYSTEM "${SLEPC_INCLUDE_DIRS}")
 else()
   if (PUGS_ENABLE_SLEPC MATCHES "^ON$")
     message(FATAL_ERROR "Could not find SLEPc!")
@@ -231,16 +225,17 @@ 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)
+  # May be risky. (make to show pugs build options)
   find_package(HDF5)
 
-  if (${HDF5_FOUND})
-    set(PUGS_CXX_FLAGS "${PUGS_CXX_FLAGS} ${HDF5_C_DEFINITIONS}")
-    include_directories(SYSTEM "${HDF5_C_INCLUDE_DIRS}")
-  else()
-    unset(HDF5_C_LIBRARIES)
-    unset(HDF5_C_HL_LIBRARIES)
-  endif()
+  # HighFive
+  set(HIGHFIVE_USE_BOOST  OFF)   # no Boost
+  set(HIGHFIVE_BUILD_DOCS OFF)   # no doc
+  set(HIGHFIVE_UNIT_TESTS OFF)   # no unit tests
+  set(HIGHFIVE_UNIT_TESTS OFF)   # no unit tests
+  set(HIGHFIVE_EXAMPLES OFF)     # no examples
+  set(HIGHFIVE_PARALLEL_HDF5 ON) # activate parallel HDF5
+  add_subdirectory(${PUGS_SOURCE_DIR}/packages/HighFive/)
 
   set(PUGS_HAS_HDF5 ${HDF5_FOUND})
 else()
@@ -306,11 +301,11 @@ endif()
 add_subdirectory("${PUGS_SOURCE_DIR}/packages/kokkos")
 
 # set as SYSTEM for static analysis
-include_directories(SYSTEM ${KOKKOS_SOURCE_DIR}/core/src)
-include_directories(SYSTEM ${KOKKOS_SOURCE_DIR}/containers/src)
-include_directories(SYSTEM ${KOKKOS_SOURCE_DIR}/tpls/desul/include)
-include_directories(SYSTEM ${KOKKOS_BINARY_DIR}/core/src)
-include_directories(SYSTEM ${KOKKOS_BINARY_DIR})
+include_directories(SYSTEM "${KOKKOS_SOURCE_DIR}/core/src")
+include_directories(SYSTEM "${KOKKOS_SOURCE_DIR}/containers/src")
+include_directories(SYSTEM "${KOKKOS_SOURCE_DIR}/tpls/desul/include")
+include_directories(SYSTEM "${KOKKOS_BINARY_DIR}/core/src")
+include_directories(SYSTEM "${KOKKOS_BINARY_DIR}")
 
 set(PUGS_BUILD_KOKKOS_DEVICES "")
 if(${Kokkos_ENABLE_PTHREAD})
@@ -362,28 +357,18 @@ endif()
 #------------------------------------------------------
 
 # Rang (colors? Useless thus necessary!)
-include_directories(${PUGS_SOURCE_DIR}/packages/rang/include)
+include_directories(SYSTEM "${PUGS_SOURCE_DIR}/packages/rang/include")
 
 # CLI11
-include_directories(${PUGS_SOURCE_DIR}/packages/CLI11/include)
+include_directories(SYSTEM "${PUGS_SOURCE_DIR}/packages/CLI11/include")
 
 # PEGTL
-include_directories(SYSTEM ${PUGS_SOURCE_DIR}/packages/PEGTL/include/tao)
-
-
-# HighFive
-set(HIGHFIVE_USE_BOOST  OFF)   # no Boost
-set(HIGHFIVE_BUILD_DOCS OFF)   # no doc
-set(HIGHFIVE_UNIT_TESTS OFF)   # no unit tests
-set(HIGHFIVE_UNIT_TESTS OFF)   # no unit tests
-set(HIGHFIVE_EXAMPLES OFF)     # no examples
-set(HIGHFIVE_PARALLEL_HDF5 ON) # activate parallel HDF5
-add_subdirectory(${PUGS_SOURCE_DIR}/packages/HighFive)
+include_directories(SYSTEM "${PUGS_SOURCE_DIR}/packages/PEGTL/include/tao")
 
 # Pugs src
-add_subdirectory(${PUGS_SOURCE_DIR}/src)
-include_directories(${PUGS_SOURCE_DIR}/src)
-include_directories(${PUGS_BINARY_DIR}/src)
+add_subdirectory("${PUGS_SOURCE_DIR}/src")
+include_directories("${PUGS_SOURCE_DIR}/src")
+include_directories("${PUGS_BINARY_DIR}/src")
 
 # Pugs tests
 set(CATCH_MODULE_PATH "${PUGS_SOURCE_DIR}/packages/Catch2")
@@ -610,6 +595,7 @@ target_link_libraries(
   PugsMesh
   PugsAlgebra
   PugsAnalysis
+  PugsDevUtils
   PugsUtils
   PugsLanguage
   PugsLanguageAST
@@ -629,8 +615,6 @@ target_link_libraries(
   ${KOKKOS_CXX_FLAGS}
   ${OPENMP_LINK_FLAGS}
   ${PUGS_STD_LINK_FLAGS}
-  ${HDF5_C_LIBRARIES}
-  ${HDF5_C_HL_LIBRARIES}
   HighFive
   stdc++fs
   )
@@ -649,6 +633,7 @@ install(TARGETS
   PugsMesh
   PugsAlgebra
   PugsAnalysis
+  PugsDevUtils
   PugsUtils
   PugsLanguage
   PugsLanguageAST
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 54cae625b..1e642e86a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -6,6 +6,9 @@ include_directories(${CMAKE_CURRENT_BINARY_DIR})
 # Pugs utils
 add_subdirectory(utils)
 
+# Pugs dev_utils
+add_subdirectory(dev_utils)
+
 # Pugs language
 add_subdirectory(language)
 
diff --git a/src/dev_utils/CMakeLists.txt b/src/dev_utils/CMakeLists.txt
new file mode 100644
index 000000000..7d709b9f7
--- /dev/null
+++ b/src/dev_utils/CMakeLists.txt
@@ -0,0 +1,11 @@
+# ------------------- Source files --------------------
+
+add_library(
+  PugsDevUtils
+  ParallelChecker.cpp
+)
+
+target_link_libraries(
+  PugsDevUtils
+  HighFive
+)
diff --git a/src/dev_utils/ParallelChecker.cpp b/src/dev_utils/ParallelChecker.cpp
new file mode 100644
index 000000000..104d30c14
--- /dev/null
+++ b/src/dev_utils/ParallelChecker.cpp
@@ -0,0 +1,35 @@
+#include <dev_utils/ParallelChecker.hpp>
+
+ParallelChecker* ParallelChecker::m_instance = nullptr;
+
+void
+ParallelChecker::create()
+{
+  Assert(ParallelChecker::m_instance == nullptr, "ParallelChecker has already been created");
+  ParallelChecker::m_instance = new ParallelChecker;
+}
+
+void
+ParallelChecker::destroy()
+{
+  Assert(ParallelChecker::m_instance != nullptr, "ParallelChecker has already been destroyed");
+  delete ParallelChecker::m_instance;
+}
+
+void
+parallel_check(const ItemValueVariant& item_value_variant,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  std::visit([&](auto&& item_value) { parallel_check(item_value, name, source_location); },
+             item_value_variant.itemValue());
+}
+
+void
+parallel_check(const DiscreteFunctionVariant& discrete_function_variant,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  std::visit([&](auto&& discrete_function) { parallel_check(discrete_function, name, source_location); },
+             discrete_function_variant.discreteFunction());
+}
diff --git a/src/dev_utils/ParallelChecker.hpp b/src/dev_utils/ParallelChecker.hpp
new file mode 100644
index 000000000..22253202a
--- /dev/null
+++ b/src/dev_utils/ParallelChecker.hpp
@@ -0,0 +1,741 @@
+#ifndef PARALLEL_CHECKER_HPP
+#define PARALLEL_CHECKER_HPP
+
+#include <highfive/H5File.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemArrayVariant.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <utils/Demangle.hpp>
+#include <utils/Filesystem.hpp>
+#include <utils/Messenger.hpp>
+#include <utils/SourceLocation.hpp>
+
+#include <fstream>
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+void parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+void parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+#ifdef PUGS_HAS_HDF5
+
+class ParallelChecker
+{
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  friend void parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
+                             const std::string& name,
+                             const SourceLocation& source_location);
+
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  friend void parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
+                             const std::string& name,
+                             const SourceLocation& source_location);
+
+ private:
+  template <typename T>
+  struct TinyVectorDataType;
+
+  template <size_t Dimension, typename DataT>
+  struct TinyVectorDataType<TinyVector<Dimension, DataT>> : public HighFive::DataType
+  {
+    TinyVectorDataType()
+    {
+      hsize_t dim[]     = {Dimension};
+      auto h5_data_type = HighFive::create_datatype<DataT>();
+      _hid              = H5Tarray_create(h5_data_type.getId(), 1, dim);
+    }
+  };
+
+  template <typename T>
+  struct TinyMatrixDataType;
+
+  template <size_t M, size_t N, typename DataT>
+  struct TinyMatrixDataType<TinyMatrix<M, N, DataT>> : public HighFive::DataType
+  {
+    TinyMatrixDataType()
+    {
+      hsize_t dim[]     = {M, N};
+      auto h5_data_type = HighFive::create_datatype<DataT>();
+      _hid              = H5Tarray_create(h5_data_type.getId(), 2, dim);
+    }
+  };
+
+  static ParallelChecker* m_instance;
+
+  size_t m_tag = 0;
+
+  std::string m_filename = "testme/parallel_checker.h5";
+
+  ParallelChecker() = default;
+
+  std::unique_ptr<HighFive::SilenceHDF5> m_silence_hdf5 = std::make_unique<HighFive::SilenceHDF5>(true);
+
+  HighFive::File
+  _createOrOpenFileRW() const
+  {
+    if (m_tag == 0) {
+      createDirectoryIfNeeded(m_filename);
+      return HighFive::File{m_filename, HighFive::File::Truncate};
+    } else {
+      return HighFive::File{m_filename, HighFive::File::ReadWrite};
+    }
+  }
+
+  void
+  _printHeader(const std::string& name, const SourceLocation& source_location) const
+  {
+    std::cout << rang::fg::cyan << " - " << rang::fgB::cyan << "parallel checker" << rang::fg::cyan << " for \""
+              << rang::fgB::magenta << name << rang::fg::cyan << "\" tag " << rang::fgB::blue << m_tag
+              << rang::fg::reset << '\n';
+    std::cout << rang::fg::cyan << " | from " << rang::fgB::blue << source_location.filename() << rang::fg::reset << ':'
+              << rang::style::bold << source_location.line() << rang::style::reset << '\n';
+  }
+
+  template <typename DataType>
+  void
+  _writeArray(HighFive::Group& group, const std::string& name, const Array<DataType>& array) const
+  {
+    using data_type = std::remove_const_t<DataType>;
+    if constexpr (is_tiny_vector_v<data_type>) {
+      auto dataset = group.createDataSet(name, HighFive::DataSpace{std::vector<size_t>{array.size()}},
+                                         TinyVectorDataType<data_type>{});
+
+      dataset.template write_raw<typename data_type::data_type>(&(array[0][0]), TinyVectorDataType<data_type>{});
+    } else if constexpr (is_tiny_matrix_v<data_type>) {
+      auto dataset = group.createDataSet(name, HighFive::DataSpace{std::vector<size_t>{array.size()}},
+                                         TinyMatrixDataType<data_type>{});
+
+      dataset.template write_raw<typename data_type::data_type>(&(array[0](0, 0)), TinyMatrixDataType<data_type>{});
+    } else {
+      auto dataset = group.createDataSet<data_type>(name, HighFive::DataSpace{std::vector<size_t>{array.size()}});
+      dataset.template write_raw<data_type>(&(array[0]));
+    }
+  }
+
+  template <typename DataType>
+  void
+  _writeTable(HighFive::Group& group, const std::string& name, const Table<DataType>& table) const
+  {
+    using data_type = std::remove_const_t<DataType>;
+    if constexpr (is_tiny_vector_v<data_type>) {
+      auto dataset =
+        group.createDataSet(name,
+                            HighFive::DataSpace{std::vector<size_t>{table.numberOfRows(), table.numberOfColumns()}},
+                            TinyVectorDataType<data_type>{});
+
+      dataset.template write_raw<typename data_type::data_type>(&(table(0, 0)[0]), TinyVectorDataType<data_type>{});
+    } else if constexpr (is_tiny_matrix_v<data_type>) {
+      auto dataset =
+        group.createDataSet(name,
+                            HighFive::DataSpace{std::vector<size_t>{table.numberOfRows(), table.numberOfColumns()}},
+                            TinyMatrixDataType<data_type>{});
+
+      dataset.template write_raw<typename data_type::data_type>(&(table(0, 0)(0, 0)), TinyMatrixDataType<data_type>{});
+    } else {
+      auto dataset =
+        group.createDataSet<data_type>(name, HighFive::DataSpace{
+                                               std::vector<size_t>{table.numberOfRows(), table.numberOfColumns()}});
+      dataset.template write_raw<data_type>(&(table(0, 0)));
+    }
+  }
+
+  template <typename DataType>
+  Array<std::remove_const_t<DataType>>
+  _readArray(HighFive::Group& group, const std::string& name) const
+  {
+    using data_type = std::remove_const_t<DataType>;
+
+    auto dataset = group.getDataSet(name);
+    Array<data_type> array(dataset.getElementCount());
+
+    if constexpr (is_tiny_vector_v<data_type>) {
+      dataset.read<data_type>(&(array[0]), TinyVectorDataType<data_type>{});
+    } else if constexpr (is_tiny_matrix_v<data_type>) {
+      dataset.read<data_type>(&(array[0]), TinyMatrixDataType<data_type>{});
+    } else {
+      dataset.read<data_type>(&(array[0]));
+    }
+    return array;
+  }
+
+  template <typename DataType>
+  Table<std::remove_const_t<DataType>>
+  _readTable(HighFive::Group& group, const std::string& name) const
+  {
+    using data_type = std::remove_const_t<DataType>;
+
+    auto dataset = group.getDataSet(name);
+    Table<data_type> table(dataset.getDimensions()[0], dataset.getDimensions()[1]);
+
+    if constexpr (is_tiny_vector_v<data_type>) {
+      dataset.read<data_type>(&(table(0, 0)), TinyVectorDataType<data_type>{});
+    } else if constexpr (is_tiny_matrix_v<data_type>) {
+      dataset.read<data_type>(&(table(0, 0)), TinyMatrixDataType<data_type>{});
+    } else {
+      dataset.read<data_type>(&(table(0, 0)));
+    }
+    return table;
+  }
+
+  template <ItemType item_type>
+  Array<const int>
+  _getItemNumber(const std::shared_ptr<const IConnectivity>& i_connectivity) const
+  {
+    switch (i_connectivity->dimension()) {
+    case 1: {
+      const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
+      return connectivity.number<item_type>().arrayView();
+    }
+    case 2: {
+      const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
+      return connectivity.number<item_type>().arrayView();
+    }
+    case 3: {
+      const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
+      return connectivity.number<item_type>().arrayView();
+    }
+    default: {
+      throw UnexpectedError("unexpected connectivity dimension");
+    }
+    }
+  }
+
+  template <ItemType item_type>
+  Array<const int>
+  _getItemOwner(const std::shared_ptr<const IConnectivity>& i_connectivity) const
+  {
+    switch (i_connectivity->dimension()) {
+    case 1: {
+      const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
+      return connectivity.owner<item_type>().arrayView();
+    }
+    case 2: {
+      const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
+      return connectivity.owner<item_type>().arrayView();
+    }
+    case 3: {
+      const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
+      return connectivity.owner<item_type>().arrayView();
+    }
+    default: {
+      throw UnexpectedError("unexpected connectivity dimension");
+    }
+    }
+  }
+
+  template <ItemType item_type>
+  void
+  _writeItemNumbers(const std::shared_ptr<const IConnectivity> i_connectivity, HighFive::Group group) const
+  {
+    this->_writeArray(group, "numbers", this->_getItemNumber<item_type>(i_connectivity));
+  }
+
+  template <typename DataType, ItemType item_type>
+  void
+  _checkIsComparable(const std::string& name,
+                     const SourceLocation& source_location,
+                     const std::vector<size_t> data_shape,
+                     const std::shared_ptr<const IConnectivity>& i_connectivity,
+                     HighFive::Group group) const
+  {
+    const std::string reference_name          = group.getAttribute("name").read<std::string>();
+    const std::string reference_file_name     = group.getAttribute("filename").read<std::string>();
+    const std::string reference_function_name = group.getAttribute("function").read<std::string>();
+    const size_t reference_line_number        = group.getAttribute("line").read<size_t>();
+    const size_t reference_dimension          = group.getAttribute("dimension").read<size_t>();
+    const std::string reference_item_type     = group.getAttribute("item_type").read<std::string>();
+    const std::string reference_data_type     = group.getAttribute("data_type").read<std::string>();
+
+    bool is_comparable = true;
+    if (i_connectivity->dimension() != reference_dimension) {
+      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different support dimensions: reference ("
+                << rang::fgB::yellow << reference_dimension << rang::fgB::red << ") / target (" << rang::fgB::yellow
+                << i_connectivity->dimension() << rang::fg::reset << ")\n";
+      is_comparable = false;
+    }
+    if (itemName(item_type) != reference_item_type) {
+      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different item types: reference (" << rang::fgB::yellow
+                << reference_item_type << rang::fgB::red << ") / target (" << rang::fgB::yellow << itemName(item_type)
+                << rang::fg::reset << ")\n";
+      is_comparable = false;
+    }
+    if (demangle<DataType>() != reference_data_type) {
+      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data types: reference (" << rang::fgB::yellow
+                << reference_data_type << rang::fgB::red << ") / target (" << rang::fgB::yellow << demangle<DataType>()
+                << rang::fg::reset << ")\n";
+      is_comparable = false;
+    }
+    std::vector reference_data_shape = group.getDataSet(reference_name).getSpace().getDimensions();
+    if (reference_data_shape.size() != data_shape.size()) {
+      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data shape kind: reference ("
+                << rang::fgB::yellow << reference_data_shape.size() << "d array" << rang::fgB::red << ") / target ("
+                << rang::fgB::yellow << data_shape.size() << "d array" << rang::fg::reset << ")\n";
+      is_comparable = false;
+    }
+    {
+      bool same_shapes = true;
+      for (size_t i = 1; i < reference_data_shape.size(); ++i) {
+        same_shapes &= (reference_data_shape[i] == data_shape[i]);
+      }
+      if (not same_shapes) {
+        std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data shape: reference ("
+                  << rang::fgB::yellow << "*";
+        for (size_t i = 1; i < reference_data_shape.size(); ++i) {
+          std::cout << ":" << reference_data_shape[i];
+        }
+        std::cout << rang::fgB::red << ") / target (" << rang::fgB::yellow << "*";
+        for (size_t i = 1; i < reference_data_shape.size(); ++i) {
+          std::cout << ":" << data_shape[i];
+        }
+        std::cout << rang::fg::reset << ")\n";
+        is_comparable = false;
+      }
+    }
+    if (name != reference_name) {
+      // Just warn for different labels (maybe useful for some kind of
+      // debugging...)
+      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "different names: reference (" << rang::fgB::yellow
+                << reference_name << rang::fgB::magenta << ") / target (" << rang::fgB::yellow << name
+                << rang::fg::reset << ")\n";
+      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "reference from " << rang::fgB::blue
+                << reference_file_name << rang::fg::reset << ':' << rang::style::bold << reference_line_number
+                << rang::style::reset << '\n';
+      if ((reference_function_name.size() > 0) or (source_location.function().size() > 0)) {
+        std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "reference function " << rang::fgB::blue
+                  << reference_function_name << rang::fg::reset << '\n';
+        std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "target function " << rang::fgB::blue
+                  << source_location.function() << rang::fg::reset << '\n';
+      }
+    }
+
+    if (not parallel::allReduceAnd(is_comparable)) {
+      throw NormalError("cannot compare data");
+    }
+  }
+
+ public:
+  static void create();
+  static void destroy();
+
+  static ParallelChecker&
+  instance()
+  {
+    return *m_instance;
+  }
+
+ private:
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  void
+  write(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
+        const std::string& name,
+        const SourceLocation& source_location)
+  {
+    this->_printHeader(name, source_location);
+
+    try {
+      HighFive::File file = this->_createOrOpenFileRW();
+
+      auto group = file.createGroup("values/" + std::to_string(m_tag));
+
+      group.createAttribute("filename", std::string{source_location.filename()});
+      group.createAttribute("function", std::string{source_location.function()});
+      group.createAttribute("line", static_cast<size_t>(source_location.line()));
+      group.createAttribute("name", name);
+
+      std::shared_ptr<const IConnectivity> i_connectivity = item_value.connectivity_ptr();
+      group.createAttribute("dimension", static_cast<size_t>(i_connectivity->dimension()));
+      group.createAttribute("item_type", std::string{itemName(item_type)});
+      group.createAttribute("data_type", demangle<DataType>());
+
+      this->_writeArray(group, name, item_value.arrayView());
+
+      this->_writeItemNumbers<item_type>(i_connectivity, group);
+
+      ++m_tag;
+
+      std::cout << rang::fg::cyan << " | writing " << rang::fgB::green << "success" << rang::fg::reset << '\n';
+    }
+    catch (HighFive::Exception& e) {
+      throw NormalError(e.what());
+    }
+  }
+
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  void
+  write(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
+        const std::string& name,
+        const SourceLocation& source_location)
+  {
+    this->_printHeader(name, source_location);
+
+    try {
+      HighFive::File file = this->_createOrOpenFileRW();
+
+      auto group = file.createGroup("values/" + std::to_string(m_tag));
+
+      group.createAttribute("filename", std::string{source_location.filename()});
+      group.createAttribute("function", std::string{source_location.function()});
+      group.createAttribute("line", static_cast<size_t>(source_location.line()));
+      group.createAttribute("name", name);
+
+      std::shared_ptr<const IConnectivity> i_connectivity = item_array.connectivity_ptr();
+      group.createAttribute("dimension", static_cast<size_t>(i_connectivity->dimension()));
+      group.createAttribute("item_type", std::string{itemName(item_type)});
+      group.createAttribute("data_type", demangle<DataType>());
+
+      this->_writeTable(group, name, item_array.tableView());
+
+      this->_writeItemNumbers<item_type>(i_connectivity, group);
+
+      ++m_tag;
+
+      std::cout << rang::fg::cyan << " | writing " << rang::fgB::green << "success" << rang::fg::reset << '\n';
+    }
+    catch (HighFive::Exception& e) {
+      throw NormalError(e.what());
+    }
+  }
+
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  void
+  compare(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
+          const std::string& name,
+          const SourceLocation& source_location)
+  {
+    this->_printHeader(name, source_location);
+
+    try {
+      HighFive::File file{m_filename, HighFive::File::ReadOnly};
+
+      auto group = file.getGroup("values/" + std::to_string(m_tag));
+
+      const std::string reference_name = group.getAttribute("name").read<std::string>();
+
+      std::shared_ptr<const IConnectivity> i_connectivity = item_value.connectivity_ptr();
+
+      this->_checkIsComparable<DataType, item_type>(name, source_location,
+                                                    std::vector<size_t>{item_value.numberOfItems()}, i_connectivity,
+                                                    group);
+
+      Array<const int> reference_item_numbers = this->_readArray<int>(group, "numbers");
+
+      Array<const DataType> reference_item_value = this->_readArray<DataType>(group, reference_name);
+
+      Array<const int> item_numbers = this->_getItemNumber<item_type>(i_connectivity);
+
+      using ItemId = ItemIdT<item_type>;
+
+      std::unordered_map<int, ItemId> item_number_to_item_id_map;
+
+      for (ItemId item_id = 0; item_id < item_numbers.size(); ++item_id) {
+        const auto& [iterator, success] =
+          item_number_to_item_id_map.insert(std::make_pair(item_numbers[item_id], item_id));
+
+        if (not success) {
+          throw UnexpectedError("item numbers have duplicate values");
+        }
+      }
+
+      Assert(item_number_to_item_id_map.size() == item_numbers.size());
+
+      Array<int> index_in_reference(item_numbers.size());
+      index_in_reference.fill(-1);
+      for (size_t i = 0; i < reference_item_numbers.size(); ++i) {
+        const auto& i_number_to_item_id = item_number_to_item_id_map.find(reference_item_numbers[i]);
+        if (i_number_to_item_id != item_number_to_item_id_map.end()) {
+          index_in_reference[i_number_to_item_id->second] = i;
+        }
+      }
+
+      if (parallel::allReduceMin(min(index_in_reference)) < 0) {
+        throw NormalError("some item numbers are not defined in reference");
+      }
+
+      Array<const int> owner = this->_getItemOwner<item_type>(i_connectivity);
+
+      bool has_own_differences = false;
+      bool is_same             = true;
+
+      for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
+        if (reference_item_value[index_in_reference[item_id]] != item_value[item_id]) {
+          is_same = false;
+          if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
+            has_own_differences = true;
+          }
+        }
+      }
+
+      is_same             = parallel::allReduceAnd(is_same);
+      has_own_differences = parallel::allReduceOr(has_own_differences);
+
+      if (is_same) {
+        std::cout << rang::fg::cyan << " | compare: " << rang::fgB::green << "success" << rang::fg::reset << '\n';
+      } else {
+        if (has_own_differences) {
+          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::red << "failed!" << rang::fg::reset;
+        } else {
+          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::yellow << "not synchronized" << rang::fg::reset;
+        }
+        std::cout << rang::fg::cyan << " [see \"" << rang::fgB::blue << "parallel_differences_" << m_tag << "_*"
+                  << rang::fg::cyan << "\" files for details]" << rang::fg::reset << '\n';
+
+        {
+          std::ofstream fout(std::string{"parallel_differences_"} + stringify(m_tag) + std::string{"_"} +
+                             stringify(parallel::rank()));
+
+          fout.precision(15);
+          for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
+            if (reference_item_value[index_in_reference[item_id]] != item_value[item_id]) {
+              const bool is_own_difference = (parallel::rank() == static_cast<size_t>(owner[item_id]));
+              if (is_own_difference) {
+                fout << rang::fgB::red << "[ own ]" << rang::fg::reset;
+              } else {
+                fout << rang::fgB::yellow << "[ghost]" << rang::fg::reset;
+              }
+              fout << " rank=" << parallel::rank() << " owner=" << owner[item_id] << " item_id=" << item_id
+                   << " number=" << item_numbers[item_id]
+                   << " reference=" << reference_item_value[index_in_reference[item_id]]
+                   << " target=" << item_value[item_id]
+                   << " difference=" << reference_item_value[index_in_reference[item_id]] - item_value[item_id] << '\n';
+              if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
+                has_own_differences = true;
+              }
+            }
+          }
+        }
+
+        if (parallel::allReduceAnd(has_own_differences)) {
+          throw NormalError("calculations differ!");
+        }
+      }
+
+      ++m_tag;
+    }
+    catch (HighFive::Exception& e) {
+      throw NormalError(e.what());
+    }
+  }
+
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  void
+  compare(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
+          const std::string& name,
+          const SourceLocation& source_location)
+  {
+    this->_printHeader(name, source_location);
+
+    try {
+      HighFive::File file{m_filename, HighFive::File::ReadOnly};
+
+      auto group = file.getGroup("values/" + std::to_string(m_tag));
+
+      const std::string reference_name = group.getAttribute("name").read<std::string>();
+
+      std::shared_ptr<const IConnectivity> i_connectivity = item_array.connectivity_ptr();
+
+      this->_checkIsComparable<DataType, item_type>(name, source_location,
+                                                    std::vector<size_t>{item_array.numberOfItems(),
+                                                                        item_array.sizeOfArrays()},
+                                                    i_connectivity, group);
+
+      Array<const int> reference_item_numbers    = this->_readArray<int>(group, "numbers");
+      Table<const DataType> reference_item_array = this->_readTable<DataType>(group, reference_name);
+
+      Array<const int> item_numbers = this->_getItemNumber<item_type>(i_connectivity);
+
+      using ItemId = ItemIdT<item_type>;
+
+      std::unordered_map<int, ItemId> item_number_to_item_id_map;
+
+      for (ItemId item_id = 0; item_id < item_numbers.size(); ++item_id) {
+        const auto& [iterator, success] =
+          item_number_to_item_id_map.insert(std::make_pair(item_numbers[item_id], item_id));
+
+        if (not success) {
+          throw UnexpectedError("item numbers have duplicate values");
+        }
+      }
+
+      Assert(item_number_to_item_id_map.size() == item_numbers.size());
+
+      Array<int> index_in_reference(item_numbers.size());
+      index_in_reference.fill(-1);
+      for (size_t i = 0; i < reference_item_numbers.size(); ++i) {
+        const auto& i_number_to_item_id = item_number_to_item_id_map.find(reference_item_numbers[i]);
+        if (i_number_to_item_id != item_number_to_item_id_map.end()) {
+          index_in_reference[i_number_to_item_id->second] = i;
+        }
+      }
+
+      if (parallel::allReduceMin(min(index_in_reference)) < 0) {
+        throw NormalError("some item numbers are not defined in reference");
+      }
+
+      Array<const int> owner = this->_getItemOwner<item_type>(i_connectivity);
+
+      bool has_own_differences = false;
+      bool is_same             = true;
+
+      for (ItemId item_id = 0; item_id < item_array.numberOfItems(); ++item_id) {
+        for (size_t i = 0; i < reference_item_array.numberOfColumns(); ++i) {
+          if (reference_item_array[index_in_reference[item_id]][i] != item_array[item_id][i]) {
+            is_same = false;
+            if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
+              has_own_differences = true;
+            }
+          }
+        }
+      }
+
+      is_same             = parallel::allReduceAnd(is_same);
+      has_own_differences = parallel::allReduceOr(has_own_differences);
+
+      if (is_same) {
+        std::cout << rang::fg::cyan << " | compare: " << rang::fgB::green << "success" << rang::fg::reset << '\n';
+      } else {
+        if (has_own_differences) {
+          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::red << "failed!" << rang::fg::reset;
+        } else {
+          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::yellow << "not synchronized" << rang::fg::reset;
+        }
+        std::cout << rang::fg::cyan << " [see \"" << rang::fgB::blue << "parallel_differences_" << m_tag << "_*"
+                  << rang::fg::cyan << "\" files for details]" << rang::fg::reset << '\n';
+
+        {
+          std::ofstream fout(std::string{"parallel_differences_"} + stringify(m_tag) + std::string{"_"} +
+                             stringify(parallel::rank()));
+
+          fout.precision(15);
+          for (ItemId item_id = 0; item_id < item_array.numberOfItems(); ++item_id) {
+            for (size_t i = 0; i < item_array.sizeOfArrays(); ++i) {
+              if (reference_item_array[index_in_reference[item_id]][i] != item_array[item_id][i]) {
+                const bool is_own_difference = (parallel::rank() == static_cast<size_t>(owner[item_id]));
+                if (is_own_difference) {
+                  fout << rang::fgB::red << "[ own ]" << rang::fg::reset;
+                } else {
+                  fout << rang::fgB::yellow << "[ghost]" << rang::fg::reset;
+                }
+                fout << " rank=" << parallel::rank() << " owner=" << owner[item_id] << " item_id=" << item_id
+                     << " column=" << i << " number=" << item_numbers[item_id]
+                     << " reference=" << reference_item_array[index_in_reference[item_id]][i]
+                     << " target=" << item_array[item_id][i]
+                     << " difference=" << reference_item_array[index_in_reference[item_id]][i] - item_array[item_id][i]
+                     << '\n';
+                if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
+                  has_own_differences = true;
+                }
+              }
+            }
+          }
+        }
+
+        if (parallel::allReduceAnd(has_own_differences)) {
+          throw NormalError("calculations differ!");
+        }
+      }
+
+      ++m_tag;
+    }
+    catch (HighFive::Exception& e) {
+      throw NormalError(e.what());
+    }
+  }
+};
+
+#else   // PUGS_HAS_HDF5
+
+class ParallelChecker
+{
+ private:
+  static ParallelChecker* m_instance;
+
+ public:
+  template <typename T>
+  void
+  write(const T&, const std::string&, const SourceLocation&)
+  {
+    throw UnexpectedError("parallel checker cannot be used without HDF5 support");
+  }
+
+  template <typename T>
+  void
+  compare(const T&, const std::string&, const SourceLocation&)
+  {
+    throw UnexpectedError("parallel checker cannot be used without HDF5 support");
+  }
+
+  static void create();
+  static void destroy();
+
+  static ParallelChecker&
+  instance()
+  {
+    return *m_instance;
+  }
+};
+
+#endif   // PUGS_HAS_HDF5
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+void
+parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  const bool write_mode = (parallel::size() == 1);
+
+  if (write_mode) {
+    ParallelChecker::instance().write(item_array, name, source_location);
+  } else {
+    ParallelChecker::instance().compare(item_array, name, source_location);
+  }
+}
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+void
+parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  const bool write_mode = (parallel::size() == 1);
+
+  if (write_mode) {
+    ParallelChecker::instance().write(item_value, name, source_location);
+  } else {
+    ParallelChecker::instance().compare(item_value, name, source_location);
+  }
+}
+
+template <size_t Dimension, typename DataType>
+void PUGS_INLINE
+parallel_check(const DiscreteFunctionP0<Dimension, DataType>& discrete_function,
+               const std::string& name,
+               const SourceLocation& source_location = SourceLocation{})
+{
+  parallel_check(discrete_function.cellValues(), name, source_location);
+}
+
+template <size_t Dimension, typename DataType>
+void PUGS_INLINE
+parallel_check(const DiscreteFunctionP0Vector<Dimension, DataType>& discrete_function,
+               const std::string& name,
+               const SourceLocation& source_location = SourceLocation{})
+{
+  parallel_check(discrete_function.cellArrays(), name, source_location);
+}
+
+void parallel_check(const ItemValueVariant& item_value_variant,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+void parallel_check(const DiscreteFunctionVariant& discrete_function_variant,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+#endif   // PARALLEL_CHECKER_HPP
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 47f490580..2e92ff898 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -1,6 +1,7 @@
 # ------------------- Source files --------------------
 
-add_library(PugsLanguageModules
+add_library(
+  PugsLanguageModules
   BinaryOperatorRegisterForVh.cpp
   BuiltinModule.cpp
   CoreModule.cpp
@@ -16,8 +17,14 @@ add_library(PugsLanguageModules
   WriterModule.cpp
 )
 
+target_link_libraries(
+  PugsLanguageModules
+  HighFive
+)
 
-add_dependencies(PugsLanguageModules
+add_dependencies(
+  PugsLanguageModules
   PugsLanguageAlgorithms
   PugsUtils
-  PugsMesh)
+  PugsMesh
+)
diff --git a/src/language/modules/DevUtilsModule.cpp b/src/language/modules/DevUtilsModule.cpp
index 2dace9adf..3d8ae43d7 100644
--- a/src/language/modules/DevUtilsModule.cpp
+++ b/src/language/modules/DevUtilsModule.cpp
@@ -1,5 +1,6 @@
 #include <language/modules/DevUtilsModule.hpp>
 
+#include <dev_utils/ParallelChecker.hpp>
 #include <language/utils/ASTDotPrinter.hpp>
 #include <language/utils/ASTExecutionInfo.hpp>
 #include <language/utils/ASTPrinter.hpp>
@@ -8,6 +9,21 @@
 
 #include <fstream>
 
+class DiscreteFunctionVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const DiscreteFunctionVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh");
+
+class ItemValueVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemValueVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_value");
+
+class ItemArrayVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemArrayVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_array");
+
 DevUtilsModule::DevUtilsModule()
 {
   this->_addBuiltinFunction("getAST", std::function(
@@ -63,6 +79,25 @@ DevUtilsModule::DevUtilsModule()
                                                 }
 
                                                 ));
+
+  this->_addBuiltinFunction("parallel_check",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function,
+                                 const std::string& name) {
+                                parallel_check(*discrete_function, name, ASTBacktrace::getInstance().sourceLocation());
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("parallel_check",
+                            std::function(
+
+                              [](const std::shared_ptr<const ItemValueVariant>& item_value, const std::string& name) {
+                                parallel_check(*item_value, name, ASTBacktrace::getInstance().sourceLocation());
+                              }
+
+                              ));
 }
 
 void
diff --git a/src/language/utils/CMakeLists.txt b/src/language/utils/CMakeLists.txt
index db3e7856b..9cf233850 100644
--- a/src/language/utils/CMakeLists.txt
+++ b/src/language/utils/CMakeLists.txt
@@ -38,10 +38,15 @@ add_library(PugsLanguageUtils
   UnaryOperatorRegisterForRn.cpp
   UnaryOperatorRegisterForRnxn.cpp
   UnaryOperatorRegisterForZ.cpp
-  )
+)
 
-
-
-add_dependencies(PugsLanguageModules
+add_dependencies(PugsLanguageUtils
+  PugsLanguageModules
   PugsUtils
-  PugsMesh)
+  PugsMesh
+)
+
+target_link_libraries(
+  PugsLanguageUtils
+  HighFive
+)
diff --git a/src/main.cpp b/src/main.cpp
index 05f01f290..b89a36ba7 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,9 +1,9 @@
 #include <analysis/QuadratureManager.hpp>
+#include <dev_utils/ParallelChecker.hpp>
 #include <language/PugsParser.hpp>
 #include <mesh/DualConnectivityManager.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
-#include <mesh/ParallelChecker.hpp>
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/PugsUtils.hpp>
 #include <utils/RandomEngine.hpp>
@@ -19,11 +19,11 @@ main(int argc, char* argv[])
   MeshDataManager::create();
   DualConnectivityManager::create();
   DualMeshManager::create();
-  parallel::ParallelChecker::create();
+  ParallelChecker::create();
 
   parser(filename);
 
-  parallel::ParallelChecker::destroy();
+  ParallelChecker::destroy();
   DualMeshManager::destroy();
   DualConnectivityManager::destroy();
   MeshDataManager::destroy();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 84c102d33..55dffcca0 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -40,5 +40,10 @@ add_library(
   MeshRandomizer.cpp
   MeshSmoother.cpp
   MeshTransformer.cpp
-  ParallelChecker.cpp
-  SynchronizerManager.cpp)
+  SynchronizerManager.cpp
+)
+
+target_link_libraries(
+  PugsMesh
+  HighFive
+)
diff --git a/src/mesh/ParallelChecker.cpp b/src/mesh/ParallelChecker.cpp
deleted file mode 100644
index 4d329e1ca..000000000
--- a/src/mesh/ParallelChecker.cpp
+++ /dev/null
@@ -1,21 +0,0 @@
-#include <mesh/ParallelChecker.hpp>
-
-namespace parallel
-{
-ParallelChecker* ParallelChecker::m_instance = nullptr;
-
-void
-ParallelChecker::create()
-{
-  Assert(ParallelChecker::m_instance == nullptr, "ParallelChecker has already been created");
-  ParallelChecker::m_instance = new ParallelChecker;
-}
-
-void
-ParallelChecker::destroy()
-{
-  Assert(ParallelChecker::m_instance != nullptr, "ParallelChecker has already been destroyed");
-  delete ParallelChecker::m_instance;
-}
-
-}   // namespace parallel
diff --git a/src/mesh/ParallelChecker.hpp b/src/mesh/ParallelChecker.hpp
deleted file mode 100644
index 8f40766bf..000000000
--- a/src/mesh/ParallelChecker.hpp
+++ /dev/null
@@ -1,363 +0,0 @@
-#ifndef PARALLEL_CHECKER_HPP
-#define PARALLEL_CHECKER_HPP
-
-#include <mesh/Connectivity.hpp>
-#include <mesh/ItemValue.hpp>
-#include <utils/HDF5.hpp>
-#include <utils/Messenger.hpp>
-
-#include <utils/SourceLocation.hpp>
-
-#include <fstream>
-#include <utils/Demangle.hpp>
-
-namespace parallel
-{
-#ifdef PUGS_HAS_HDF5
-
-template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-void check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
-           const std::string& name,
-           const SourceLocation& source_location = SourceLocation{});
-
-class ParallelChecker
-{
- private:
-  static ParallelChecker* m_instance;
-
-  size_t m_tag = 0;
-
-  std::string m_filename = "parallel_checker.h5";
-
-  ParallelChecker() = default;
-
-  void
-  _printHeader(const std::string& name, const SourceLocation& source_location) const
-  {
-    std::cout << rang::fg::cyan << " | " << rang::fgB::cyan << "parallel checker" << rang::fg::cyan << " for \""
-              << rang::fgB::magenta << name << rang::fg::cyan << "\" tag " << rang::fgB::blue << m_tag
-              << rang::fg::reset << '\n';
-    std::cout << rang::fg::cyan << " | from " << rang::fgB::blue << source_location.filename() << rang::fg::reset << ':'
-              << rang::style::bold << source_location.line() << rang::style::reset << '\n';
-  }
-
- public:
-  static void create();
-  static void destroy();
-
-  static ParallelChecker&
-  instance()
-  {
-    return *m_instance;
-  }
-
-  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-  friend void check(const ItemValue<DataType, item_type, ConnectivityPtr>&, const std::string&, const SourceLocation&);
-
- private:
-  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-  void
-  write(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
-        const std::string& name,
-        const SourceLocation& source_location)
-  {
-    this->_printHeader(name, source_location);
-
-    auto file_id = [&] {
-      if (m_tag == 0) {
-        return HDF5::create(m_filename);
-      } else {
-        return HDF5::openFileRW(m_filename);
-      }
-    }();
-
-    auto values_group_id = HDF5::createOrOpenGroup(file_id, "/values");
-    auto group_id        = HDF5::createOrOpenGroup(values_group_id, std::to_string(m_tag));
-
-    HDF5::writeAttribute(group_id, "filename", std::string{source_location.filename()});
-    HDF5::writeAttribute(group_id, "function", source_location.function());
-    HDF5::writeAttribute(group_id, "line", static_cast<size_t>(source_location.line()));
-    HDF5::writeAttribute(group_id, "name", name);
-
-    std::shared_ptr<const IConnectivity> i_connectivity = item_value.connectivity_ptr();
-    HDF5::writeAttribute(group_id, "dimension", static_cast<size_t>(i_connectivity->dimension()));
-    HDF5::writeAttribute(group_id, "item_type", itemName(item_type));
-    HDF5::writeAttribute(group_id, "data_type", demangle<DataType>());
-
-    HDF5::writeArray(group_id, name, item_value.arrayView());
-
-    switch (i_connectivity->dimension()) {
-    case 1: {
-      const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
-      HDF5::writeArray(group_id, "numbers", connectivity.number<item_type>().arrayView());
-      break;
-    }
-    case 2: {
-      const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
-      HDF5::writeArray(group_id, "numbers", connectivity.number<item_type>().arrayView());
-      break;
-    }
-    case 3: {
-      const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
-      HDF5::writeArray(group_id, "numbers", connectivity.number<item_type>().arrayView());
-      break;
-    }
-    default: {
-      throw UnexpectedError("unexpected connectivity dimension");
-    }
-    }
-
-    ++m_tag;
-
-    HDF5::close(values_group_id);
-    HDF5::close(group_id);
-    HDF5::close(file_id);
-
-    std::cout << rang::fg::cyan << " | writing " << rang::fgB::green << "success" << rang::fg::reset << '\n';
-  }
-
-  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-  void
-  compare(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
-          const std::string& name,
-          const SourceLocation& source_location)
-  {
-    this->_printHeader(name, source_location);
-
-    auto file_id = HDF5::openFileRO(m_filename);
-
-    auto values_group_id = HDF5::openGroup(file_id, "/values");
-    auto group_id        = HDF5::openGroup(values_group_id, std::to_string(m_tag));
-
-    const std::string reference_name          = HDF5::readAttribute<std::string>(group_id, "name");
-    const std::string reference_file_name     = HDF5::readAttribute<std::string>(group_id, "filename");
-    const std::string reference_function_name = HDF5::readAttribute<std::string>(group_id, "function");
-    const size_t reference_line_number        = HDF5::readAttribute<size_t>(group_id, "line");
-    const size_t reference_dimension          = HDF5::readAttribute<size_t>(group_id, "dimension");
-    const std::string reference_item_type     = HDF5::readAttribute<std::string>(group_id, "item_type");
-    const std::string reference_data_type     = HDF5::readAttribute<std::string>(group_id, "data_type");
-
-    std::shared_ptr<const IConnectivity> i_connectivity = item_value.connectivity_ptr();
-
-    bool is_comparable = true;
-    if (i_connectivity->dimension() != reference_dimension) {
-      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different support dimensions: reference ("
-                << rang::fgB::yellow << reference_dimension << rang::fgB::red << ") / target (" << rang::fgB::yellow
-                << i_connectivity->dimension() << rang::fg::reset << ")\n";
-      is_comparable = false;
-    }
-    if (itemName(item_type) != reference_item_type) {
-      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different item types: reference (" << rang::fgB::yellow
-                << reference_item_type << rang::fgB::red << ") / target (" << rang::fgB::yellow << itemName(item_type)
-                << rang::fg::reset << ")\n";
-      is_comparable = false;
-    }
-    if (demangle<DataType>() != reference_data_type) {
-      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data types: reference (" << rang::fgB::yellow
-                << reference_data_type << rang::fgB::red << ") / target (" << rang::fgB::yellow << demangle<DataType>()
-                << rang::fg::reset << ")\n";
-      is_comparable = false;
-    }
-    if (name != reference_name) {
-      // Just warn for different labels (maybe useful for some kind of
-      // debugging...)
-      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "different names: reference (" << rang::fgB::yellow
-                << reference_name << rang::fgB::magenta << ") / target (" << rang::fgB::yellow << name
-                << rang::fg::reset << ")\n";
-      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "reference from " << rang::fgB::blue
-                << reference_file_name << rang::fg::reset << ':' << rang::style::bold << reference_line_number
-                << rang::style::reset << '\n';
-      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "reference function " << rang::fgB::blue
-                << reference_function_name << rang::fg::reset << '\n';
-      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "target function " << rang::fgB::blue
-                << source_location.function() << rang::fg::reset << '\n';
-    }
-
-    if (not parallel::allReduceAnd(is_comparable)) {
-      throw NormalError("cannot compare data");
-    }
-
-    Array<const int> reference_item_numbers = HDF5::readArray<int>(group_id, "numbers");
-    Array<const DataType> reference_item_value =
-      HDF5::readArray<std::remove_const_t<DataType> >(group_id, reference_name);
-
-    Array<const int> item_numbers = [&] {
-      switch (i_connectivity->dimension()) {
-      case 1: {
-        const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
-        return connectivity.number<item_type>().arrayView();
-      }
-      case 2: {
-        const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
-        return connectivity.number<item_type>().arrayView();
-      }
-      case 3: {
-        const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
-        return connectivity.number<item_type>().arrayView();
-      }
-      default: {
-        throw UnexpectedError("unexpected connectivity dimension");
-      }
-      }
-    }();
-
-    using ItemId = ItemIdT<item_type>;
-
-    std::unordered_map<int, ItemId> item_number_to_item_id_map;
-
-    for (ItemId item_id = 0; item_id < item_numbers.size(); ++item_id) {
-      const auto& [iterator, success] =
-        item_number_to_item_id_map.insert(std::make_pair(item_numbers[item_id], item_id));
-
-      if (not success) {
-        throw UnexpectedError("item numbers have duplicate values");
-      }
-    }
-
-    Assert(item_number_to_item_id_map.size() == item_numbers.size());
-
-    Array<int> index_in_reference(item_numbers.size());
-    index_in_reference.fill(-1);
-    for (size_t i = 0; i < reference_item_numbers.size(); ++i) {
-      const auto& i_number_to_item_id = item_number_to_item_id_map.find(reference_item_numbers[i]);
-      if (i_number_to_item_id != item_number_to_item_id_map.end()) {
-        index_in_reference[i_number_to_item_id->second] = i;
-      }
-    }
-
-    if (parallel::allReduceMin(min(index_in_reference)) < 0) {
-      throw NormalError("some item numbers are not defined in reference");
-    }
-
-    Array<const int> owner = [&] {
-      switch (i_connectivity->dimension()) {
-      case 1: {
-        const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
-        return connectivity.owner<item_type>().arrayView();
-      }
-      case 2: {
-        const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
-        return connectivity.owner<item_type>().arrayView();
-      }
-      case 3: {
-        const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
-        return connectivity.owner<item_type>().arrayView();
-      }
-      default: {
-        throw UnexpectedError("unexpected connectivity dimension");
-      }
-      }
-    }();
-
-    bool has_own_differences = false;
-    bool is_same             = true;
-
-    for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
-      if (reference_item_value[index_in_reference[item_id]] != item_value[item_id]) {
-        is_same = false;
-        if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
-          has_own_differences = true;
-        }
-      }
-    }
-
-    is_same             = parallel::allReduceAnd(is_same);
-    has_own_differences = parallel::allReduceOr(has_own_differences);
-
-    if (is_same) {
-      std::cout << rang::fg::cyan << " | compare: " << rang::fgB::green << "success" << rang::fg::reset << '\n';
-    } else {
-      if (has_own_differences) {
-        std::cout << rang::fg::cyan << " | compare: " << rang::fgB::red << "failed!" << rang::fg::reset;
-      } else {
-        std::cout << rang::fg::cyan << " | compare: " << rang::fgB::yellow << "not synchronized" << rang::fg::reset;
-      }
-      std::cout << rang::fg::cyan << " [see \"" << rang::fgB::blue << "parallel_differences_" << m_tag << "_*"
-                << rang::fg::cyan << "\" files for details]" << rang::fg::reset << '\n';
-
-      {
-        std::ofstream fout(std::string{"parallel_differences_"} + stringify(m_tag) + std::string{"_"} +
-                           stringify(parallel::rank()));
-
-        fout.precision(15);
-        for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
-          if (reference_item_value[index_in_reference[item_id]] != item_value[item_id]) {
-            const bool is_own_difference = (parallel::rank() == static_cast<size_t>(owner[item_id]));
-            if (is_own_difference) {
-              fout << rang::fgB::red << "[ own ]" << rang::fg::reset;
-            } else {
-              fout << rang::fgB::yellow << "[ghost]" << rang::fg::reset;
-            }
-            fout << " rank=" << parallel::rank() << " owner=" << owner[item_id] << " item_id=" << item_id
-                 << " number=" << item_numbers[item_id]
-                 << " reference=" << reference_item_value[index_in_reference[item_id]]
-                 << " target=" << item_value[item_id]
-                 << " difference=" << reference_item_value[index_in_reference[item_id]] - item_value[item_id] << '\n';
-            if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
-              has_own_differences = true;
-            }
-          }
-        }
-      }
-
-      if (parallel::allReduceAnd(has_own_differences)) {
-        throw NormalError("calculations differ!");
-      }
-    }
-
-    HDF5::close(values_group_id);
-    HDF5::close(group_id);
-    HDF5::close(file_id);
-    ++m_tag;
-  }
-};
-
-template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-void
-check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
-      const std::string& name,
-      const SourceLocation& source_location)
-{
-  const bool write_mode = (parallel::size() == 1);
-
-  std::cout << '\n';
-  if (write_mode) {
-    ParallelChecker::instance().write(item_value, name, source_location);
-  } else {
-    ParallelChecker::instance().compare(item_value, name, source_location);
-  }
-  std::cout << '\n';
-}
-
-#else   // PUGS_HAS_HDF5
-
-template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-void
-check(const ItemValue<DataType, item_type, ConnectivityPtr>&,
-      const std::string&,
-      const SourceLocation& = SourceLocation{})
-{
-  throw UnexpectedError("parallel checker cannot be used without HDF5 support");
-}
-
-class ParallelChecker
-{
- private:
-  static ParallelChecker* m_instance;
-
- public:
-  static void create();
-  static void destroy();
-
-  static ParallelChecker&
-  instance()
-  {
-    return *m_instance;
-  }
-};
-
-#endif   // PUGS_HAS_HDF5
-
-}   // namespace parallel
-
-#endif   // PARALLEL_CHECKER_HPP
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index f3173dd26..12b22d01f 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -6,3 +6,8 @@ add_library(
   GnuplotWriter1D.cpp
   VTKWriter.cpp
   WriterBase.cpp)
+
+target_link_libraries(
+  PugsOutput
+  HighFive
+)
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 1379aea67..5d49a9177 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -9,4 +9,10 @@ add_library(
   DiscreteFunctionUtils.cpp
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
-  FluxingAdvectionSolver.cpp)
+  FluxingAdvectionSolver.cpp
+)
+
+target_link_libraries(
+  PugsScheme
+  HighFive
+)
diff --git a/src/scheme/DiscreteFunctionVariant.hpp b/src/scheme/DiscreteFunctionVariant.hpp
index 6feb86bcb..af1b1cc37 100644
--- a/src/scheme/DiscreteFunctionVariant.hpp
+++ b/src/scheme/DiscreteFunctionVariant.hpp
@@ -92,7 +92,7 @@ class DiscreteFunctionVariant
                   "DiscreteFunctionP0Vector with this DataType is not allowed in variant");
   }
 
-  DiscreteFunctionVariant& operator=(DiscreteFunctionVariant&&) = default;
+  DiscreteFunctionVariant& operator=(DiscreteFunctionVariant&&)      = default;
   DiscreteFunctionVariant& operator=(const DiscreteFunctionVariant&) = default;
 
   DiscreteFunctionVariant(const DiscreteFunctionVariant&) = default;
diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 5c922089e..05b002e7c 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -18,7 +18,7 @@
 #endif   // PUGS_HAS_PETSC
 
 #ifdef PUGS_HAS_HDF5
-#include <hdf5.h>
+#include <highfive/H5File.hpp>
 #endif   // PUGS_HAS_HDF5
 
 std::string
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 9284d3f03..8a32abb36 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -9,7 +9,6 @@ add_library(
   Demangle.cpp
   Exceptions.cpp
   FPEManager.cpp
-  HDF5.cpp
   Messenger.cpp
   Partitioner.cpp
   PETScWrapper.cpp
@@ -25,6 +24,7 @@ target_link_libraries(
   PugsUtils
   ${PETSC_LIBRARIES}
   ${SLEPC_LIBRARIES}
+  HighFive
 )
 
 # --------------- get git revision info ---------------
diff --git a/src/utils/HDF5.cpp b/src/utils/HDF5.cpp
deleted file mode 100644
index 1361ddbe9..000000000
--- a/src/utils/HDF5.cpp
+++ /dev/null
@@ -1,235 +0,0 @@
-#include <utils/HDF5.hpp>
-
-#ifdef PUGS_HAS_HDF5
-
-#include <utils/Exceptions.hpp>
-
-#include <fstream>
-#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 < 0) {
-    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 < 0) {
-    std::ostringstream error_msg;
-    error_msg << "cannot open file '" << rang::fgB::yellow << filename << rang::fg::reset << "' for reading";
-    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::createOrOpenGroup(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 open group '" << rang::fgB::yellow << groupname << rang::fg::reset << "'";
-    throw NormalError(error_msg.str());
-  }
-}
-
-HDF5::GroupId
-HDF5::openGroup(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 open group '" << rang::fgB::yellow << groupname << rang::fg::reset << "'";
-    throw NormalError(error_msg.str());
-  }
-}
-
-template <HDF5::IdType id_type>
-void
-HDF5::_writeArray(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");   // LCOV_EXCL_LINE
-  }
-  herr_t status = H5Dwrite(dataset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, array);
-  if (status < 0) {
-    throw NormalError("cannot write dataset");   // LCOV_EXCL_LINE
-  }
-
-  H5Dclose(dataset);
-  H5Sclose(dataspace);
-}
-
-template void HDF5::_writeArray(const FileId& id,
-                                const std::string& name,
-                                const DatatypeId& type,
-                                const size_t& array_size,
-                                const void* array);
-
-template void HDF5::_writeArray(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::_writeTable(const HDFId<id_type>& id,
-                  const std::string& name,
-                  const DatatypeId& type,
-                  const size_t& table_nb_rows,
-                  const size_t& table_nb_columns,
-                  const void* table)
-{
-  hsize_t dataspace_size[] = {table_nb_rows, table_nb_columns};
-  hid_t dataspace          = H5Screate_simple(2, 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");   // LCOV_EXCL_LINE
-  }
-  herr_t status = H5Dwrite(dataset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, table);
-  if (status < 0) {
-    throw NormalError("cannot write dataset");   // LCOV_EXCL_LINE
-  }
-
-  H5Dclose(dataset);
-  H5Sclose(dataspace);
-}
-
-template void HDF5::_writeTable(const FileId& id,
-                                const std::string& name,
-                                const DatatypeId& type,
-                                const size_t& table_nb_rows,
-                                const size_t& table_nb_columns,
-                                const void* table);
-
-template void HDF5::_writeTable(const GroupId& id,
-                                const std::string& name,
-                                const DatatypeId& type,
-                                const size_t& table_nb_rows,
-                                const size_t& table_nb_columns,
-                                const void* table);
-
-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);
-
-  if (H5Aexists(id, name.c_str()) > 0) {
-    std::ostringstream error_msg;
-    error_msg << "attribute '" << rang::fgB::yellow << name << rang::fg::reset << "' already exists";
-    throw NormalError(error_msg.str());
-  }
-
-  hid_t attribute_id = H5Acreate2(id, name.c_str(), type, space, H5P_DEFAULT, H5P_DEFAULT);
-  if (attribute_id < 0) {
-  }
-
-  H5Awrite(attribute_id, 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 << "'";
-      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
deleted file mode 100644
index 933d1658d..000000000
--- a/src/utils/HDF5.hpp
+++ /dev/null
@@ -1,436 +0,0 @@
-#ifndef HDF5_HPP
-#define HDF5_HPP
-
-#include <utils/pugs_config.hpp>
-
-#ifdef PUGS_HAS_HDF5
-
-#include <utils/Array.hpp>
-#include <utils/Exceptions.hpp>
-#include <utils/Table.hpp>
-
-#include <hdf5.h>
-
-#include <utils/Demangle.hpp>
-
-class HDF5
-{
- private:
-  enum class IdType
-  {
-    attribute,
-    dataset,
-    file,
-    group,
-    datatype
-  };
-
-  template <IdType type>
-  class [[nodiscard]] 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 _writeArray(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 _writeTable(const HDFId<id_type>& id,
-                          const std::string& data_name,
-                          const DatatypeId& data_type,
-                          const size_t& table_nb_rows,
-                          const size_t& table_nb_columns,
-                          const void* table);
-
-  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 if constexpr (is_tiny_vector_v<std::decay_t<T>>) {
-      using ContentT      = typename T::data_type;
-      hsize_t dim[]       = {DataT::Dimension};
-      DatatypeId datatype = H5Tarray_create(_getDataType<ContentT>(), 1, dim);
-      return datatype;
-    } else if constexpr (is_tiny_matrix_v<std::decay_t<T>>) {
-      using ContentT      = typename T::data_type;
-      hsize_t dim[]       = {DataT::NumberOfRows, DataT::NumberOfColumns};
-      DatatypeId datatype = H5Tarray_create(_getDataType<ContentT>(), 2, dim);
-      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 createOrOpenGroup(const GroupId& group_id, const std::string& group_name);
-
-  static GroupId openGroup(const FileId& file_id, const std::string& group_name);
-  static GroupId openGroup(const GroupId& group_id, const std::string& group_name);
-
-  template <typename DataT, IdType id_type>
-  static void
-  writeArray(const HDFId<id_type>& id, const std::string& name, const Array<DataT>& array)
-  {
-    if (H5Lexists(id, name.c_str(), H5P_DEFAULT)) {
-      std::ostringstream error_msg;
-      error_msg << "dataset '" << name << "' already exists";
-      throw NormalError(error_msg.str());
-    }
-
-    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)]);
-
-    _writeArray(id, name, datatype, array.size(), array_start);
-
-    close(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);
-
-    if (not H5Lexists(id, name.c_str(), H5P_DEFAULT)) {
-      std::ostringstream error_msg;
-      error_msg << "cannot open dataset '" << name << "'";
-      throw NormalError(error_msg.str());
-    }
-
-    DatasetId dataset_id = H5Dopen2(id, name.c_str(), H5P_DEFAULT);
-
-    hid_t dataset_type = H5Dget_type(dataset_id);
-    if (not H5Tequal(dataset_type, datatype)) {
-      std::ostringstream error_msg;
-      error_msg << "invalid type for array '" << rang::fgB::yellow << name << rang::fg::reset << "'";
-      throw NormalError(error_msg.str());
-    }
-
-    hid_t dataspace_id = H5Dget_space(dataset_id);
-    int ndims          = H5Sget_simple_extent_ndims(dataspace_id);
-    if (ndims != 1) {
-      throw NormalError("stored dataset is not an array");
-    }
-    Array<hsize_t> dims(ndims);
-    H5Sget_simple_extent_dims(dataspace_id, &(dims[0]), NULL);
-
-    Array<DataT> array{dims[0]};
-
-    H5Dread(dataset_id, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(array[0]));
-
-    close(dataset_id);
-    close(datatype);
-
-    return array;
-  }
-
-  template <typename DataT, IdType id_type>
-  static void
-  writeTable(const HDFId<id_type>& id, const std::string& name, const Table<DataT>& table)
-  {
-    if (H5Lexists(id, name.c_str(), H5P_DEFAULT)) {
-      std::ostringstream error_msg;
-      error_msg << "dataset '" << name << "' already exists";
-      throw NormalError(error_msg.str());
-    }
-
-    static_assert((id_type == IdType::file) or (id_type == IdType::group));
-    DatatypeId datatype = _getDataType<DataT>();
-    H5Tset_order(datatype, H5T_ORDER_LE);
-
-    const void* table_start = &(table(typename Table<DataT>::index_type(0), typename Table<DataT>::index_type(0)));
-
-    _writeTable(id, name, datatype, table.numberOfRows(), table.numberOfColumns(), table_start);
-
-    close(datatype);
-  }
-
-  template <typename DataT, IdType id_type>
-  static Table<DataT>
-  readTable(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);
-
-    if (not H5Lexists(id, name.c_str(), H5P_DEFAULT)) {
-      std::ostringstream error_msg;
-      error_msg << "cannot open dataset '" << name << "'";
-      throw NormalError(error_msg.str());
-    }
-
-    DatasetId dataset_id = H5Dopen2(id, name.c_str(), H5P_DEFAULT);
-
-    hid_t dataset_type = H5Dget_type(dataset_id);
-    if (not H5Tequal(dataset_type, datatype)) {
-      std::ostringstream error_msg;
-      error_msg << "invalid type for table '" << rang::fgB::yellow << name << rang::fg::reset << "'";
-      throw NormalError(error_msg.str());
-    }
-
-    hid_t dataspace_id = H5Dget_space(dataset_id);
-    int ndims          = H5Sget_simple_extent_ndims(dataspace_id);
-    if (ndims != 2) {
-      throw NormalError("stored dataset is not a table");
-    }
-    Array<hsize_t> dims(ndims);
-    H5Sget_simple_extent_dims(dataspace_id, &(dims[0]), NULL);
-
-    Table<DataT> table{dims[0], dims[1]};
-
-    H5Dread(dataset_id, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(table(0, 0)));
-
-    close(dataset_id);
-    close(datatype);
-
-    return table;
-  }
-
-  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);
-      close(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);
-      close(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);
-      close(datatype);
-    } else if constexpr (std::is_arithmetic_v<std::decay_t<DataT>> or is_tiny_vector_v<std::decay_t<DataT>> or
-                         is_tiny_matrix_v<std::decay_t<DataT>>) {
-      DatatypeId datatype = _getDataType<DataT>();
-      _writeAttribute(id, name, datatype, &value);
-      close(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);
-
-      close(datatype);
-      return value;
-    } else if constexpr (std::is_arithmetic_v<std::decay_t<DataT>> or is_tiny_vector_v<DataT> or
-                         is_tiny_matrix_v<DataT>) {
-      DatatypeId datatype = _getDataType<DataT>();
-      DataT value;
-      _readAttribute(id, name, datatype, &value);
-      close(datatype);
-      return value;
-    } else {
-      static_assert(std::is_same_v<std::decay_t<DataT>, std::string>, "unexpected attribute type");
-    }
-  }
-
-  static void
-  close(const AttributeId& attribute_id)
-  {
-    H5Aclose(attribute_id);
-  }
-
-  static void
-  close(const DatasetId& dataset_id)
-  {
-    H5Dclose(dataset_id);
-  }
-
-  static void
-  close(const DatatypeId& file_id)
-  {
-    H5Tclose(file_id);
-  }
-
-  static void
-  close(const FileId& file_id)
-  {
-    H5Fclose(file_id);
-  }
-
-  static void
-  close(const GroupId& file_id)
-  {
-    H5Gclose(file_id);
-  }
-};
-#endif   // PUGS_HAS_HDF5
-
-#endif   // HDF5_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 6985929d0..c47e88837 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -96,7 +96,6 @@ 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
@@ -226,14 +225,14 @@ target_link_libraries (unit_tests
   PugsScheme
   PugsOutput
   PugsUtils
+  PugsDevUtils
   Kokkos::kokkos
   ${PARMETIS_LIBRARIES}
   ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
   ${PETSC_LIBRARIES}
   Catch2
   ${PUGS_STD_LINK_FLAGS}
-  ${HDF5_C_LIBRARIES}
-  ${HDF5_C_HL_LIBRARIES}
+  HighFive
   stdc++fs
   )
 
@@ -252,6 +251,7 @@ target_link_libraries (mpi_unit_tests
   PugsLanguageUtils
   PugsScheme
   PugsOutput
+  PugsDevUtils
   PugsUtils
   PugsAlgebra
   PugsMesh
@@ -261,8 +261,7 @@ target_link_libraries (mpi_unit_tests
   ${PETSC_LIBRARIES}
   Catch2
   ${PUGS_STD_LINK_FLAGS}
-  ${HDF5_C_LIBRARIES}
-  ${HDF5_C_HL_LIBRARIES}
+  HighFive
   stdc++fs
   )
 
diff --git a/tests/test_HDF5.cpp b/tests/test_HDF5.cpp
deleted file mode 100644
index 2d424b063..000000000
--- a/tests/test_HDF5.cpp
+++ /dev/null
@@ -1,388 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <utils/HDF5.hpp>
-
-#include <algebra/TinyMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-
-#ifdef PUGS_HAS_HDF5
-
-#include <filesystem>
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("HDF5", "[utils]")
-{
-  H5Eset_auto(H5E_DEFAULT, NULL, NULL);
-
-  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"};
-  }();
-
-  auto is_same_array = [](auto f, auto g) -> bool {
-    if (f.size() != g.size()) {
-      return false;
-    }
-    for (size_t i = 0; i < f.size(); ++i) {
-      if (f[i] != g[i]) {
-        return false;
-      }
-    }
-    return true;
-  };
-
-  auto is_same_table = [](auto f, auto g) -> bool {
-    if ((f.numberOfRows() != g.numberOfRows()) or (f.numberOfColumns() != g.numberOfColumns())) {
-      return false;
-    }
-    for (size_t i = 0; i < f.numberOfRows(); ++i) {
-      for (size_t j = 0; j < f.numberOfColumns(); ++j) {
-        if (f(i, j) != g(i, j)) {
-          return false;
-        }
-      }
-    }
-    return true;
-  };
-
-  Array<const int> int_array = [] {
-    Array<int> array(12);
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] = 3 * i + 2;
-    }
-    return array;
-  }();
-
-  Array<const double> double_array = [] {
-    Array<double> array(10);
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] = 3.7 * i - 2.3;
-    }
-    return array;
-  }();
-
-  Array<const TinyVector<3, double>> R3_array = [] {
-    Array<TinyVector<3, double>> array(9);
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] = TinyVector<3, double>{3.7 * i - 2.3, 2.3 + i, -3.2 * i - 1};
-    }
-    return array;
-  }();
-
-  Array<const TinyMatrix<2, 3, int>> Z23_array = [] {
-    Array<TinyMatrix<2, 3, int>> array(13);
-    for (int i = 0; i < static_cast<int>(array.size()); ++i) {
-      array[i] = TinyMatrix<2, 3, int>{3 * i - 2,  2 + i,   //
-                                       -3 * i - 1, 2 + i,   //
-                                       2 * i + 3,  3 - 2 * i};
-    }
-    return array;
-  }();
-
-  Table<const int> int_table = [] {
-    Table<int> table(12, 3);
-    for (size_t i = 0; i < table.numberOfRows(); ++i) {
-      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
-        table(i, j) = 3 * i + 2 - 2 * j;
-      }
-    }
-    return table;
-  }();
-
-  Table<const double> double_table = [] {
-    Table<double> table(10, 5);
-    for (size_t i = 0; i < table.numberOfRows(); ++i) {
-      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
-        table(i, j) = 3.2 * i + 2 - 2.5 * j;
-      }
-    }
-    return table;
-  }();
-
-  SECTION("create file")
-  {
-    {
-      std::stringstream error_msg;
-      error_msg << "error: cannot open file '" << filename << "' for writing";
-      REQUIRE_THROWS_WITH(HDF5::openFileRW(filename), error_msg.str());
-    }
-
-    {
-      std::stringstream error_msg;
-      error_msg << "error: cannot open file '" << filename << "' for reading";
-      REQUIRE_THROWS_WITH(HDF5::openFileRO(filename), error_msg.str());
-    }
-
-    HDF5::FileId file_id = HDF5::create(filename);
-    REQUIRE(file_id >= 0);
-    HDF5::close(file_id);
-
-    REQUIRE_THROWS_WITH(HDF5::create("/a/non/existing/path/to/a/file"),
-                        "error: cannot create file '/a/non/existing/path/to/a/file'");
-  }
-
-  SECTION("add attributes at file root")
-  {
-    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::writeAttribute(file_id, "tinyvector attr", TinyVector<2, int>(2, 3));
-    HDF5::writeAttribute(file_id, "tinymatrix attr", TinyMatrix<4, 3, double>(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12));
-
-    REQUIRE_THROWS_WITH(HDF5::writeAttribute(file_id, "foo attr", std::string{"bar"}),
-                        "error: attribute 'foo attr' already exists");
-
-    HDF5::close(file_id);
-  }
-
-  SECTION("add data at file root")
-  {
-    HDF5::FileId file_id = HDF5::openFileRW(filename);
-    REQUIRE(file_id >= 0);
-
-    HDF5::writeArray(file_id, "int array", int_array);
-    HDF5::writeArray(file_id, "double array", double_array);
-    HDF5::writeArray(file_id, "R3 array", R3_array);
-
-    HDF5::writeTable(file_id, "int table", int_table);
-    HDF5::writeTable(file_id, "double table", double_table);
-
-    REQUIRE_THROWS_WITH(HDF5::writeArray(file_id, "double array", double_array),
-                        "error: dataset 'double array' already exists");
-
-    REQUIRE_THROWS_WITH(HDF5::writeTable(file_id, "double table", double_table),
-                        "error: dataset 'double table' already exists");
-
-    HDF5::close(file_id);
-  }
-
-  SECTION("add groups and group attributes")
-  {
-    HDF5::FileId file_id = HDF5::openFileRW(filename);
-    REQUIRE(file_id >= 0);
-
-    HDF5::GroupId group_A_id = HDF5::createOrOpenGroup(file_id, "group_A");
-    {
-      REQUIRE_NOTHROW(HDF5::createOrOpenGroup(file_id, "group_A"));
-    }
-    REQUIRE(group_A_id >= 0);
-
-    HDF5::writeAttribute(group_A_id, "char attr", static_cast<char>(-'c'));
-    HDF5::writeAttribute(group_A_id, "uchar attr", static_cast<unsigned char>('a'));
-    HDF5::writeAttribute(group_A_id, "short attr", static_cast<short>(-23));
-    HDF5::writeAttribute(group_A_id, "unsigned short attr", static_cast<unsigned short>(31));
-    HDF5::writeAttribute(group_A_id, "int attr", static_cast<int>(-121));
-    HDF5::writeAttribute(group_A_id, "unsigned int attr", static_cast<unsigned int>(263));
-    HDF5::writeAttribute(group_A_id, "long int attr", static_cast<long int>(-112));
-    HDF5::writeAttribute(group_A_id, "long unsigned int attr", static_cast<unsigned long int>(103));
-    HDF5::writeAttribute(group_A_id, "int8_t attr", static_cast<int8_t>(-215));
-    HDF5::writeAttribute(group_A_id, "uint8_t attr", static_cast<uint8_t>(237));
-    HDF5::writeAttribute(group_A_id, "int16_t attr", static_cast<int16_t>(-124));
-    HDF5::writeAttribute(group_A_id, "uint16_t attr", static_cast<uint16_t>(244));
-    HDF5::writeAttribute(group_A_id, "int32_t attr", static_cast<int32_t>(-1243));
-    HDF5::writeAttribute(group_A_id, "uint32_t attr", static_cast<uint32_t>(1427));
-    HDF5::writeAttribute(group_A_id, "int64_t attr", static_cast<int64_t>(-15156));
-    HDF5::writeAttribute(group_A_id, "uint64_t attr", static_cast<uint64_t>(19999));
-    HDF5::writeAttribute(group_A_id, "float attr", static_cast<float>(1.4142));
-    HDF5::writeAttribute(group_A_id, "double attr", static_cast<double>(1.4142135));
-    HDF5::writeAttribute(group_A_id, "long double attr", static_cast<long double>(1.41421356237));
-    HDF5::writeAttribute(group_A_id, "foo attr", std::string{"foobar"});
-
-    HDF5::writeArray(group_A_id, "int array", int_array);
-    HDF5::writeTable(group_A_id, "int table", int_table);
-
-    REQUIRE_THROWS_WITH(HDF5::writeArray(group_A_id, "int array", int_array),
-                        "error: dataset 'int array' already exists");
-    REQUIRE_THROWS_WITH(HDF5::writeTable(group_A_id, "int table", int_table),
-                        "error: dataset 'int table' already exists");
-
-    HDF5::GroupId group_B_id = HDF5::createOrOpenGroup(group_A_id, "group_B");
-    REQUIRE(group_B_id >= 0);
-    {
-      REQUIRE_NOTHROW(HDF5::createOrOpenGroup(group_A_id, "group_B"));
-    }
-
-    HDF5::writeAttribute(group_B_id, "char attr", static_cast<char>(-'e'));
-    HDF5::writeAttribute(group_B_id, "uchar attr", static_cast<unsigned char>('s'));
-    HDF5::writeAttribute(group_B_id, "short attr", static_cast<short>(-133));
-    HDF5::writeAttribute(group_B_id, "unsigned short attr", static_cast<unsigned short>(33));
-    HDF5::writeAttribute(group_B_id, "int attr", static_cast<int>(-131));
-    HDF5::writeAttribute(group_B_id, "unsigned int attr", static_cast<unsigned int>(323));
-    HDF5::writeAttribute(group_B_id, "long int attr", static_cast<long int>(-1041));
-    HDF5::writeAttribute(group_B_id, "long unsigned int attr", static_cast<unsigned long int>(1301));
-    HDF5::writeAttribute(group_B_id, "int8_t attr", static_cast<int8_t>(-325));
-    HDF5::writeAttribute(group_B_id, "uint8_t attr", static_cast<uint8_t>(267));
-    HDF5::writeAttribute(group_B_id, "int16_t attr", static_cast<int16_t>(-194));
-    HDF5::writeAttribute(group_B_id, "uint16_t attr", static_cast<uint16_t>(824));
-    HDF5::writeAttribute(group_B_id, "int32_t attr", static_cast<int32_t>(-1234));
-    HDF5::writeAttribute(group_B_id, "uint32_t attr", static_cast<uint32_t>(1437));
-    HDF5::writeAttribute(group_B_id, "int64_t attr", static_cast<int64_t>(-15415));
-    HDF5::writeAttribute(group_B_id, "uint64_t attr", static_cast<uint64_t>(18999));
-    HDF5::writeAttribute(group_B_id, "float attr", static_cast<float>(2.718));
-    HDF5::writeAttribute(group_B_id, "double attr", static_cast<double>(2.718281));
-    HDF5::writeAttribute(group_B_id, "long double attr", static_cast<long double>(2.71828182846));
-    HDF5::writeAttribute(group_B_id, "foo attr", std::string{"foobar"});
-
-    HDF5::writeArray(group_B_id, "double array", double_array);
-    HDF5::writeArray(group_B_id, "Z23 array", Z23_array);
-    REQUIRE_THROWS_WITH(HDF5::writeArray(group_B_id, "double array", double_array),
-                        "error: dataset 'double array' already exists");
-
-    HDF5::close(group_A_id);
-    HDF5::close(group_B_id);
-    HDF5::close(file_id);
-  }
-
-  SECTION("open file ro and open groups")
-  {
-    std::filesystem::permissions(filename,
-                                 std::filesystem::perms::owner_write | std::filesystem::perms::group_all |
-                                   std::filesystem::perms::others_all,
-                                 std::filesystem::perm_options::remove);
-
-    HDF5::FileId file_id = HDF5::openFileRO(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(HDF5::readAttribute<TinyVector<2, int>>(file_id, "tinyvector attr") == TinyVector<2, int>(2, 3));
-    REQUIRE(HDF5::readAttribute<TinyMatrix<4, 3, double>>(file_id, "tinymatrix attr") ==
-            TinyMatrix<4, 3, double>(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12));
-
-    REQUIRE(is_same_array(HDF5::readArray<TinyVector<3, double>>(file_id, "R3 array"), R3_array));
-    REQUIRE(is_same_array(HDF5::readArray<int>(file_id, "int array"), int_array));
-    REQUIRE(is_same_array(HDF5::readArray<double>(file_id, "double array"), double_array));
-
-    REQUIRE(is_same_table(HDF5::readTable<int>(file_id, "int table"), int_table));
-    REQUIRE(is_same_table(HDF5::readTable<double>(file_id, "double table"), double_table));
-
-    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: invalid type for attribute 'int64_t attr'");
-
-    REQUIRE_THROWS_WITH(HDF5::readArray<float>(file_id, "invalid array"), "error: cannot open dataset 'invalid array'");
-    REQUIRE_THROWS_WITH(HDF5::readArray<float>(file_id, "double array"),
-                        "error: invalid type for array 'double array'");
-
-    REQUIRE_THROWS_WITH(HDF5::readTable<float>(file_id, "invalid table"), "error: cannot open dataset 'invalid table'");
-    REQUIRE_THROWS_WITH(HDF5::readTable<float>(file_id, "double table"),
-                        "error: invalid type for table 'double table'");
-
-    REQUIRE_THROWS_WITH(HDF5::readArray<double>(file_id, "double table"), "error: stored dataset is not an array");
-    REQUIRE_THROWS_WITH(HDF5::readTable<double>(file_id, "double array"), "error: stored dataset is not a table");
-
-    HDF5::GroupId group_A_id = HDF5::openGroup(file_id, "group_A");
-    REQUIRE(group_A_id >= 0);
-
-    REQUIRE(HDF5::readAttribute<char>(group_A_id, "char attr") == static_cast<char>(-'c'));
-    REQUIRE(HDF5::readAttribute<unsigned char>(group_A_id, "uchar attr") == static_cast<unsigned char>('a'));
-    REQUIRE(HDF5::readAttribute<short>(group_A_id, "short attr") == static_cast<short>(-23));
-    REQUIRE(HDF5::readAttribute<unsigned short>(group_A_id, "unsigned short attr") == static_cast<unsigned short>(31));
-    REQUIRE(HDF5::readAttribute<int>(group_A_id, "int attr") == static_cast<int>(-121));
-    REQUIRE(HDF5::readAttribute<unsigned int>(group_A_id, "unsigned int attr") == static_cast<unsigned int>(263));
-    REQUIRE(HDF5::readAttribute<long int>(group_A_id, "long int attr") == static_cast<long int>(-112));
-    REQUIRE(HDF5::readAttribute<unsigned long int>(group_A_id, "long unsigned int attr") ==
-            static_cast<unsigned long int>(103));
-    REQUIRE(HDF5::readAttribute<int8_t>(group_A_id, "int8_t attr") == static_cast<int8_t>(-215));
-    REQUIRE(HDF5::readAttribute<uint8_t>(group_A_id, "uint8_t attr") == static_cast<uint8_t>(237));
-    REQUIRE(HDF5::readAttribute<int16_t>(group_A_id, "int16_t attr") == static_cast<int16_t>(-124));
-    REQUIRE(HDF5::readAttribute<uint16_t>(group_A_id, "uint16_t attr") == static_cast<uint16_t>(244));
-    REQUIRE(HDF5::readAttribute<int32_t>(group_A_id, "int32_t attr") == static_cast<int32_t>(-1243));
-    REQUIRE(HDF5::readAttribute<uint32_t>(group_A_id, "uint32_t attr") == static_cast<uint32_t>(1427));
-    REQUIRE(HDF5::readAttribute<int64_t>(group_A_id, "int64_t attr") == static_cast<int64_t>(-15156));
-    REQUIRE(HDF5::readAttribute<uint64_t>(group_A_id, "uint64_t attr") == static_cast<uint64_t>(19999));
-    REQUIRE(HDF5::readAttribute<float>(group_A_id, "float attr") == static_cast<float>(1.4142));
-    REQUIRE(HDF5::readAttribute<double>(group_A_id, "double attr") == static_cast<double>(1.4142135));
-    REQUIRE(HDF5::readAttribute<long double>(group_A_id, "long double attr") ==
-            static_cast<long double>(1.41421356237));
-    REQUIRE(HDF5::readAttribute<std::string>(group_A_id, "foo attr") == std::string{"foobar"});
-
-    REQUIRE(is_same_array(HDF5::readArray<int>(group_A_id, "int array"), int_array));
-    REQUIRE_THROWS_WITH(HDF5::readArray<double>(group_A_id, "double array"),
-                        "error: cannot open dataset 'double array'");
-
-    HDF5::GroupId group_B_id = HDF5::openGroup(group_A_id, "group_B");
-    REQUIRE(group_B_id >= 0);
-
-    REQUIRE_THROWS_WITH(HDF5::openGroup(file_id, "group_B"), "error: cannot open group 'group_B'");
-    REQUIRE_THROWS_WITH(HDF5::openGroup(group_A_id, "group_C"), "error: cannot open group 'group_C'");
-
-    REQUIRE(HDF5::readAttribute<char>(group_B_id, "char attr") == static_cast<char>(-'e'));
-    REQUIRE(HDF5::readAttribute<unsigned char>(group_B_id, "uchar attr") == static_cast<unsigned char>('s'));
-    REQUIRE(HDF5::readAttribute<short>(group_B_id, "short attr") == static_cast<short>(-133));
-    REQUIRE(HDF5::readAttribute<unsigned short>(group_B_id, "unsigned short attr") == static_cast<unsigned short>(33));
-    REQUIRE(HDF5::readAttribute<int>(group_B_id, "int attr") == static_cast<int>(-131));
-    REQUIRE(HDF5::readAttribute<unsigned int>(group_B_id, "unsigned int attr") == static_cast<unsigned int>(323));
-    REQUIRE(HDF5::readAttribute<long int>(group_B_id, "long int attr") == static_cast<long int>(-1041));
-    REQUIRE(HDF5::readAttribute<unsigned long int>(group_B_id, "long unsigned int attr") ==
-            static_cast<unsigned long int>(1301));
-    REQUIRE(HDF5::readAttribute<int8_t>(group_B_id, "int8_t attr") == static_cast<int8_t>(-325));
-    REQUIRE(HDF5::readAttribute<uint8_t>(group_B_id, "uint8_t attr") == static_cast<uint8_t>(267));
-    REQUIRE(HDF5::readAttribute<int16_t>(group_B_id, "int16_t attr") == static_cast<int16_t>(-194));
-    REQUIRE(HDF5::readAttribute<uint16_t>(group_B_id, "uint16_t attr") == static_cast<uint16_t>(824));
-    REQUIRE(HDF5::readAttribute<int32_t>(group_B_id, "int32_t attr") == static_cast<int32_t>(-1234));
-    REQUIRE(HDF5::readAttribute<uint32_t>(group_B_id, "uint32_t attr") == static_cast<uint32_t>(1437));
-    REQUIRE(HDF5::readAttribute<int64_t>(group_B_id, "int64_t attr") == static_cast<int64_t>(-15415));
-    REQUIRE(HDF5::readAttribute<uint64_t>(group_B_id, "uint64_t attr") == static_cast<uint64_t>(18999));
-    REQUIRE(HDF5::readAttribute<float>(group_B_id, "float attr") == static_cast<float>(2.718));
-    REQUIRE(HDF5::readAttribute<double>(group_B_id, "double attr") == static_cast<double>(2.718281));
-    REQUIRE(HDF5::readAttribute<long double>(group_B_id, "long double attr") ==
-            static_cast<long double>(2.71828182846));
-
-    REQUIRE(HDF5::readAttribute<std::string>(group_B_id, "foo attr") == std::string{"foobar"});
-
-    REQUIRE(is_same_array(HDF5::readArray<double>(group_B_id, "double array"), double_array));
-    REQUIRE(is_same_array(HDF5::readArray<TinyMatrix<2, 3, int>>(group_B_id, "Z23 array"), Z23_array));
-
-    REQUIRE_THROWS_WITH(HDF5::readArray<int>(group_B_id, "int array"), "error: cannot open dataset 'int array'");
-
-    HDF5::close(group_A_id);
-    HDF5::close(group_B_id);
-    HDF5::close(file_id);
-  }
-}
-#endif   // PUGS_HAS_HDF5
-- 
GitLab