From d36b8a0d3009bb8a9b0c3e6948baf9cea6607ad5 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 16 Oct 2023 23:03:30 +0200
Subject: [PATCH] Add parallel checker support for SubItem{Array,Value}PerItem

---
 src/dev/CMakeLists.txt                        |   5 +
 src/dev/ParallelChecker.cpp                   |  18 -
 src/dev/ParallelChecker.hpp                   | 600 +++++++++++++++++-
 ...ParallelCheckerDiscreteFunctionVariant.cpp |  10 +
 src/dev/ParallelCheckerItemArrayVariant.cpp   |  10 +
 src/dev/ParallelCheckerItemValueVariant.cpp   |  10 +
 ...allelCheckerSubItemArrayPerItemVariant.cpp |  10 +
 ...allelCheckerSubItemValuePerItemVariant.cpp |  10 +
 src/language/modules/DevUtilsModule.cpp       |  41 ++
 src/mesh/SubItemArrayPerItem.hpp              |   7 +
 src/mesh/SubItemArrayPerItemVariant.hpp       |   4 +-
 src/mesh/SubItemValuePerItem.hpp              |   8 +
 src/mesh/SubItemValuePerItemVariant.hpp       |   4 +-
 src/utils/PugsUtils.cpp                       |   4 +-
 14 files changed, 700 insertions(+), 41 deletions(-)
 create mode 100644 src/dev/ParallelCheckerDiscreteFunctionVariant.cpp
 create mode 100644 src/dev/ParallelCheckerItemArrayVariant.cpp
 create mode 100644 src/dev/ParallelCheckerItemValueVariant.cpp
 create mode 100644 src/dev/ParallelCheckerSubItemArrayPerItemVariant.cpp
 create mode 100644 src/dev/ParallelCheckerSubItemValuePerItemVariant.cpp

diff --git a/src/dev/CMakeLists.txt b/src/dev/CMakeLists.txt
index 22264f7ed..a212a76fa 100644
--- a/src/dev/CMakeLists.txt
+++ b/src/dev/CMakeLists.txt
@@ -3,6 +3,11 @@
 add_library(
   PugsDev
   ParallelChecker.cpp
+  ParallelCheckerDiscreteFunctionVariant.cpp
+  ParallelCheckerItemArrayVariant.cpp
+  ParallelCheckerItemValueVariant.cpp
+  ParallelCheckerSubItemArrayPerItemVariant.cpp
+  ParallelCheckerSubItemValuePerItemVariant.cpp
 )
 
 target_link_libraries(
diff --git a/src/dev/ParallelChecker.cpp b/src/dev/ParallelChecker.cpp
index 680cd92e6..2e6880fb8 100644
--- a/src/dev/ParallelChecker.cpp
+++ b/src/dev/ParallelChecker.cpp
@@ -15,21 +15,3 @@ 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/ParallelChecker.hpp b/src/dev/ParallelChecker.hpp
index b19689a63..e195a4280 100644
--- a/src/dev/ParallelChecker.hpp
+++ b/src/dev/ParallelChecker.hpp
@@ -9,6 +9,8 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemValueVariant.hpp>
+#include <mesh/SubItemArrayPerItemVariant.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Demangle.hpp>
 #include <utils/Filesystem.hpp>
@@ -23,7 +25,17 @@ void parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_
                     const SourceLocation& source_location = SourceLocation{});
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
-void parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value,
+void parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+void parallel_check(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_value_per_item,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+void parallel_check(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_array_per_item,
                     const std::string& name,
                     const SourceLocation& source_location = SourceLocation{});
 
@@ -58,6 +70,16 @@ class ParallelChecker
                              const std::string& name,
                              const SourceLocation& source_location);
 
+  template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  friend void parallel_check(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_value_per_item,
+                             const std::string& name,
+                             const SourceLocation& source_location);
+
+  template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  friend void parallel_check(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_array_per_item,
+                             const std::string& name,
+                             const SourceLocation& source_location);
+
 #ifdef PUGS_HAS_HDF5
  private:
   template <typename T>
@@ -237,6 +259,29 @@ class ParallelChecker
     }
   }
 
+  template <ItemType item_type, ItemType sub_item_type>
+  Array<const typename ConnectivityMatrix::IndexType>
+  _getSubItemRowsMap(const std::shared_ptr<const IConnectivity>& i_connectivity)
+  {
+    switch (i_connectivity->dimension()) {
+    case 1: {
+      const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
+      return connectivity.getMatrix(item_type, sub_item_type).rowsMap();
+    }
+    case 2: {
+      const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
+      return connectivity.getMatrix(item_type, sub_item_type).rowsMap();
+    }
+    case 3: {
+      const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
+      return connectivity.getMatrix(item_type, sub_item_type).rowsMap();
+    }
+    default: {
+      throw UnexpectedError("unexpected connectivity dimension");
+    }
+    }
+  }
+
   template <ItemType item_type>
   Array<const int>
   _getItemOwner(const std::shared_ptr<const IConnectivity>& i_connectivity) const
@@ -278,8 +323,31 @@ class ParallelChecker
     group.createHardLink("numbers", item_numbers);
   }
 
-  template <typename DataType, ItemType item_type>
+  template <typename ItemOfItem>
   void
+  _writeSubItemRowsMap(const std::shared_ptr<const IConnectivity> i_connectivity,
+                       HighFive::File file,
+                       HighFive::Group group)
+  {
+    constexpr ItemType item_type     = ItemOfItem::item_type;
+    constexpr ItemType sub_item_type = ItemOfItem::sub_item_type;
+
+    std::string subitem_of_item_row_map_path =
+      "/connectivities/" + std::to_string(this->_getConnectivityId(i_connectivity)) + '/' +
+      std::string{itemName(sub_item_type)} + "of" + std::string{itemName(item_type)};
+
+    if (not file.exist(subitem_of_item_row_map_path)) {
+      HighFive::Group subitem_of_item_row_map_group = file.createGroup(subitem_of_item_row_map_path);
+      this->_writeArray(subitem_of_item_row_map_group, "rows_map",
+                        this->_getSubItemRowsMap<item_type, sub_item_type>(i_connectivity));
+    }
+
+    HighFive::DataSet subitem_of_item_row_map = file.getDataSet(subitem_of_item_row_map_path + "/rows_map");
+    group.createHardLink("rows_map", subitem_of_item_row_map);
+  }
+
+  template <typename DataType, ItemType item_type>
+  bool
   _checkIsComparable(const std::string& name,
                      const SourceLocation& source_location,
                      const std::vector<size_t> data_shape,
@@ -343,8 +411,8 @@ class ParallelChecker
       // 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";
+                << rang::style::reversed << reference_name << rang::style::reset << rang::fgB::magenta << ") / target ("
+                << rang::fgB::yellow << rang::style::reversed << name << rang::style::reset << 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';
@@ -356,6 +424,53 @@ class ParallelChecker
       }
     }
 
+    return is_comparable;
+  }
+
+  template <ItemType sub_item_type>
+  bool
+  _checkIsComparable(HighFive::Group group) const
+  {
+    const std::string reference_sub_item_type = group.getAttribute("subitem_type").read<std::string>();
+
+    bool is_comparable = true;
+    if (itemName(sub_item_type) != reference_sub_item_type) {
+      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different sub_item types: reference ("
+                << rang::fgB::yellow << reference_sub_item_type << rang::fgB::red << ") / target (" << rang::fgB::yellow
+                << itemName(sub_item_type) << rang::fg::reset << ")\n";
+      is_comparable = false;
+    }
+
+    return is_comparable;
+  }
+
+  template <typename DataType, ItemType item_type>
+  void
+  _throwIfNotComparable(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
+  {
+    bool is_comparable =
+      this->_checkIsComparable<DataType, item_type>(name, source_location, data_shape, i_connectivity, group);
+    if (not parallel::allReduceAnd(is_comparable)) {
+      throw NormalError("cannot compare data");
+    }
+  }
+
+  template <typename DataType, typename ItemOfItem>
+  void
+  _throwIfNotComparable(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
+  {
+    bool is_comparable = this->_checkIsComparable<DataType, ItemOfItem::item_type>(name, source_location, data_shape,
+                                                                                   i_connectivity, group)   //
+                         and this->_checkIsComparable<ItemOfItem::sub_item_type>(group);
+
     if (not parallel::allReduceAnd(is_comparable)) {
       throw NormalError("cannot compare data");
     }
@@ -436,6 +551,92 @@ class ParallelChecker
     }
   }
 
+  template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  void
+  write(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_value_per_item,
+        const std::string& name,
+        const SourceLocation& source_location)
+  {
+    constexpr ItemType item_type     = ItemOfItem::item_type;
+    constexpr ItemType sub_item_type = ItemOfItem::sub_item_type;
+
+    HighFive::SilenceHDF5 m_silence_hdf5{true};
+    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 = subitem_value_per_item.connectivity_ptr();
+      group.createAttribute("dimension", static_cast<size_t>(i_connectivity->dimension()));
+      group.createAttribute("item_type", std::string{itemName(item_type)});
+      group.createAttribute("subitem_type", std::string{itemName(sub_item_type)});
+
+      group.createAttribute("data_type", demangle<DataType>());
+
+      this->_writeArray(group, name, subitem_value_per_item.arrayView());
+
+      this->_writeItemNumbers<item_type>(i_connectivity, file, group);
+      this->_writeSubItemRowsMap<ItemOfItem>(i_connectivity, file, 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, typename ItemOfItem, typename ConnectivityPtr>
+  void
+  write(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_value_per_item,
+        const std::string& name,
+        const SourceLocation& source_location)
+  {
+    constexpr ItemType item_type     = ItemOfItem::item_type;
+    constexpr ItemType sub_item_type = ItemOfItem::sub_item_type;
+
+    HighFive::SilenceHDF5 m_silence_hdf5{true};
+    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 = subitem_value_per_item.connectivity_ptr();
+      group.createAttribute("dimension", static_cast<size_t>(i_connectivity->dimension()));
+      group.createAttribute("item_type", std::string{itemName(item_type)});
+      group.createAttribute("subitem_type", std::string{itemName(sub_item_type)});
+
+      group.createAttribute("data_type", demangle<DataType>());
+
+      this->_writeTable(group, name, subitem_value_per_item.tableView());
+
+      this->_writeItemNumbers<item_type>(i_connectivity, file, group);
+      this->_writeSubItemRowsMap<ItemOfItem>(i_connectivity, file, 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,
@@ -454,9 +655,9 @@ class ParallelChecker
 
       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);
+      this->_throwIfNotComparable<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");
 
@@ -575,10 +776,10 @@ class ParallelChecker
 
       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);
+      this->_throwIfNotComparable<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);
@@ -683,20 +884,347 @@ class ParallelChecker
     }
   }
 
+  template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  void
+  compare(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_value_per_item,
+          const std::string& name,
+          const SourceLocation& source_location)
+  {
+    HighFive::SilenceHDF5 m_silence_hdf5{true};
+    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 = subitem_value_per_item.connectivity_ptr();
+
+      this->_throwIfNotComparable<DataType, ItemOfItem>(name, source_location,
+                                                        std::vector<size_t>{subitem_value_per_item.numberOfItems()},
+                                                        i_connectivity, group);
+
+      constexpr ItemType item_type     = ItemOfItem::item_type;
+      constexpr ItemType sub_item_type = ItemOfItem::sub_item_type;
+      using IndexType                  = typename ConnectivityMatrix::IndexType;
+
+      Array<const int> reference_item_numbers           = this->_readArray<int>(group, "numbers");
+      Array<const IndexType> reference_subitem_rows_map = this->_readArray<IndexType>(group, "rows_map");
+
+      Array<const DataType> reference_subitem_value_per_item = this->_readArray<DataType>(group, reference_name);
+
+      Array<const int> item_numbers           = this->_getItemNumber<item_type>(i_connectivity);
+      Array<const IndexType> sub_item_row_map = this->_getSubItemRowsMap<item_type, sub_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> item_index_in_reference(item_numbers.size());
+      item_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()) {
+          item_index_in_reference[i_number_to_item_id->second] = i;
+        }
+      }
+
+      if (parallel::allReduceMin(min(item_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 < subitem_value_per_item.numberOfItems(); ++item_id) {
+        const size_t reference_item_index     = item_index_in_reference[item_id];
+        const size_t index_begin_in_reference = reference_subitem_rows_map[reference_item_index];
+        const size_t index_end_in_reference   = reference_subitem_rows_map[reference_item_index + 1];
+
+        bool item_is_same = true;
+        if ((index_end_in_reference - index_begin_in_reference) != subitem_value_per_item[item_id].size()) {
+          item_is_same = false;
+        } else {
+          for (size_t i_sub_item = 0; i_sub_item < subitem_value_per_item[item_id].size(); ++i_sub_item) {
+            if (reference_subitem_value_per_item[index_begin_in_reference + i_sub_item] !=
+                subitem_value_per_item[item_id][i_sub_item]) {
+              item_is_same = false;
+            }
+          }
+        }
+        if (not item_is_same) {
+          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 < subitem_value_per_item.numberOfItems(); ++item_id) {
+            const size_t reference_item_index     = item_index_in_reference[item_id];
+            const size_t index_begin_in_reference = reference_subitem_rows_map[reference_item_index];
+            const size_t index_end_in_reference   = reference_subitem_rows_map[reference_item_index + 1];
+            if ((index_end_in_reference - index_begin_in_reference) != subitem_value_per_item[item_id].size()) {
+              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[subitems number]=" << index_end_in_reference - index_begin_in_reference
+                   << " target[subitems number]=" << subitem_value_per_item[item_id].size() << '\n';
+
+            } else {
+              for (size_t i_sub_item = 0; i_sub_item < subitem_value_per_item[item_id].size(); ++i_sub_item) {
+                if (reference_subitem_value_per_item[index_begin_in_reference + i_sub_item] !=
+                    subitem_value_per_item[item_id][i_sub_item]) {
+                  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] << " i_subitem=" << i_sub_item
+                       << " reference=" << reference_subitem_value_per_item[index_begin_in_reference + i_sub_item]
+                       << " target=" << subitem_value_per_item[item_id][i_sub_item] << " difference="
+                       << reference_subitem_value_per_item[index_begin_in_reference + i_sub_item] -
+                            subitem_value_per_item[item_id][i_sub_item]
+                       << '\n';
+                }
+              }
+            }
+          }
+        }
+
+        if (has_own_differences) {
+          throw NormalError("calculations differ!");
+        }
+      }
+
+      ++m_tag;
+    }
+    catch (HighFive::Exception& e) {
+      throw NormalError(e.what());
+    }
+  }
+
+  template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  void
+  compare(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_array_per_item,
+          const std::string& name,
+          const SourceLocation& source_location)
+  {
+    HighFive::SilenceHDF5 m_silence_hdf5{true};
+    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 = subitem_array_per_item.connectivity_ptr();
+
+      this->_throwIfNotComparable<DataType, ItemOfItem>(name, source_location,
+                                                        std::vector<size_t>{subitem_array_per_item.numberOfItems(),
+                                                                            subitem_array_per_item.sizeOfArrays()},
+                                                        i_connectivity, group);
+
+      constexpr ItemType item_type     = ItemOfItem::item_type;
+      constexpr ItemType sub_item_type = ItemOfItem::sub_item_type;
+      using IndexType                  = typename ConnectivityMatrix::IndexType;
+
+      Array<const int> reference_item_numbers           = this->_readArray<int>(group, "numbers");
+      Array<const IndexType> reference_subitem_rows_map = this->_readArray<IndexType>(group, "rows_map");
+
+      Table<const DataType> reference_subitem_array_per_item = this->_readTable<DataType>(group, reference_name);
+
+      Array<const int> item_numbers           = this->_getItemNumber<item_type>(i_connectivity);
+      Array<const IndexType> sub_item_row_map = this->_getSubItemRowsMap<item_type, sub_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> item_index_in_reference(item_numbers.size());
+      item_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()) {
+          item_index_in_reference[i_number_to_item_id->second] = i;
+        }
+      }
+
+      if (parallel::allReduceMin(min(item_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 < subitem_array_per_item.numberOfItems(); ++item_id) {
+        const size_t reference_item_index     = item_index_in_reference[item_id];
+        const size_t index_begin_in_reference = reference_subitem_rows_map[reference_item_index];
+        const size_t index_end_in_reference   = reference_subitem_rows_map[reference_item_index + 1];
+
+        bool item_is_same = true;
+        if ((index_end_in_reference - index_begin_in_reference) != subitem_array_per_item[item_id].numberOfRows()) {
+          item_is_same = false;
+        } else {
+          for (size_t i_sub_item = 0; i_sub_item < subitem_array_per_item[item_id].numberOfRows(); ++i_sub_item) {
+            for (size_t i = 0; i < subitem_array_per_item.sizeOfArrays(); ++i) {
+              if (reference_subitem_array_per_item[index_begin_in_reference + i_sub_item][i] !=
+                  subitem_array_per_item[item_id][i_sub_item][i]) {
+                item_is_same = false;
+              }
+            }
+          }
+        }
+        if (not item_is_same) {
+          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 < subitem_array_per_item.numberOfItems(); ++item_id) {
+            const size_t reference_item_index     = item_index_in_reference[item_id];
+            const size_t index_begin_in_reference = reference_subitem_rows_map[reference_item_index];
+            const size_t index_end_in_reference   = reference_subitem_rows_map[reference_item_index + 1];
+            if ((index_end_in_reference - index_begin_in_reference) != subitem_array_per_item[item_id].numberOfRows()) {
+              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[subitems number]=" << index_end_in_reference - index_begin_in_reference
+                   << " target[subitems number]=" << subitem_array_per_item[item_id].numberOfRows() << '\n';
+
+            } else {
+              for (size_t i_sub_item = 0; i_sub_item < subitem_array_per_item[item_id].numberOfRows(); ++i_sub_item) {
+                for (size_t i = 0; i < subitem_array_per_item.sizeOfArrays(); ++i) {
+                  if (reference_subitem_array_per_item[index_begin_in_reference + i_sub_item][i] !=
+                      subitem_array_per_item[item_id][i_sub_item][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
+                         << " number=" << item_numbers[item_id] << " i_subitem=" << i_sub_item << " i_value=" << i
+                         << " reference=" << reference_subitem_array_per_item[index_begin_in_reference + i_sub_item]
+                         << " target=" << subitem_array_per_item[item_id][i_sub_item] << " difference="
+                         << reference_subitem_array_per_item[index_begin_in_reference + i_sub_item][i] -
+                              subitem_array_per_item[item_id][i_sub_item][i]
+                         << '\n';
+                  }
+                }
+              }
+            }
+          }
+        }
+
+        if (has_own_differences) {
+          throw NormalError("calculations differ!");
+        }
+      }
+
+      ++m_tag;
+    }
+    catch (HighFive::Exception& e) {
+      throw NormalError(e.what());
+    }
+  }
+
 #else    // PUGS_HAS_HDF5
 
   template <typename T>
   void
   write(const T&, const std::string&, const SourceLocation&)
   {
-    throw UnexpectedError("parallel checker cannot be used without HDF5 support");
+    throw NormalError("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");
+    throw NormalError("parallel checker cannot be used without HDF5 support");
   }
 #endif   // PUGS_HAS_HDF5
 
@@ -767,6 +1295,19 @@ class ParallelChecker
   }
 };
 
+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)
+{
+  if (ParallelChecker::instance().isWriting()) {
+    ParallelChecker::instance().write(item_value, name, source_location);
+  } else {
+    ParallelChecker::instance().compare(item_value, name, source_location);
+  }
+}
+
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
 parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
@@ -780,16 +1321,29 @@ parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array
   }
 }
 
-template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
 void
-parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
+parallel_check(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_value_per_item,
                const std::string& name,
                const SourceLocation& source_location)
 {
   if (ParallelChecker::instance().isWriting()) {
-    ParallelChecker::instance().write(item_value, name, source_location);
+    ParallelChecker::instance().write(subitem_value_per_item, name, source_location);
   } else {
-    ParallelChecker::instance().compare(item_value, name, source_location);
+    ParallelChecker::instance().compare(subitem_value_per_item, name, source_location);
+  }
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+void
+parallel_check(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& subitem_array_per_item,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  if (ParallelChecker::instance().isWriting()) {
+    ParallelChecker::instance().write(subitem_array_per_item, name, source_location);
+  } else {
+    ParallelChecker::instance().compare(subitem_array_per_item, name, source_location);
   }
 }
 
@@ -815,6 +1369,18 @@ void parallel_check(const ItemValueVariant& item_value_variant,
                     const std::string& name,
                     const SourceLocation& source_location = SourceLocation{});
 
+void parallel_check(const ItemArrayVariant& item_array_variant,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+void parallel_check(const SubItemValuePerItemVariant& subitem_value_per_item_variant,
+                    const std::string& name,
+                    const SourceLocation& source_location = SourceLocation{});
+
+void parallel_check(const SubItemArrayPerItemVariant& subitem_array_per_item_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{});
diff --git a/src/dev/ParallelCheckerDiscreteFunctionVariant.cpp b/src/dev/ParallelCheckerDiscreteFunctionVariant.cpp
new file mode 100644
index 000000000..626fa6842
--- /dev/null
+++ b/src/dev/ParallelCheckerDiscreteFunctionVariant.cpp
@@ -0,0 +1,10 @@
+#include <dev/ParallelChecker.hpp>
+
+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/ParallelCheckerItemArrayVariant.cpp b/src/dev/ParallelCheckerItemArrayVariant.cpp
new file mode 100644
index 000000000..2ca0f7e3f
--- /dev/null
+++ b/src/dev/ParallelCheckerItemArrayVariant.cpp
@@ -0,0 +1,10 @@
+#include <dev/ParallelChecker.hpp>
+
+void
+parallel_check(const ItemArrayVariant& item_array_variant,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  std::visit([&](auto&& item_value) { parallel_check(item_value, name, source_location); },
+             item_array_variant.itemArray());
+}
diff --git a/src/dev/ParallelCheckerItemValueVariant.cpp b/src/dev/ParallelCheckerItemValueVariant.cpp
new file mode 100644
index 000000000..4cbc9ed21
--- /dev/null
+++ b/src/dev/ParallelCheckerItemValueVariant.cpp
@@ -0,0 +1,10 @@
+#include <dev/ParallelChecker.hpp>
+
+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());
+}
diff --git a/src/dev/ParallelCheckerSubItemArrayPerItemVariant.cpp b/src/dev/ParallelCheckerSubItemArrayPerItemVariant.cpp
new file mode 100644
index 000000000..d33898f50
--- /dev/null
+++ b/src/dev/ParallelCheckerSubItemArrayPerItemVariant.cpp
@@ -0,0 +1,10 @@
+#include <dev/ParallelChecker.hpp>
+
+void
+parallel_check(const SubItemArrayPerItemVariant& subitem_array_per_item_variant,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  std::visit([&](auto&& item_value) { parallel_check(item_value, name, source_location); },
+             subitem_array_per_item_variant.subItemArrayPerItem());
+}
diff --git a/src/dev/ParallelCheckerSubItemValuePerItemVariant.cpp b/src/dev/ParallelCheckerSubItemValuePerItemVariant.cpp
new file mode 100644
index 000000000..4584202e2
--- /dev/null
+++ b/src/dev/ParallelCheckerSubItemValuePerItemVariant.cpp
@@ -0,0 +1,10 @@
+#include <dev/ParallelChecker.hpp>
+
+void
+parallel_check(const SubItemValuePerItemVariant& subitem_value_per_item_variant,
+               const std::string& name,
+               const SourceLocation& source_location)
+{
+  std::visit([&](auto&& item_value) { parallel_check(item_value, name, source_location); },
+             subitem_value_per_item_variant.subItemValuePerItem());
+}
diff --git a/src/language/modules/DevUtilsModule.cpp b/src/language/modules/DevUtilsModule.cpp
index 1de94fa63..7685443b0 100644
--- a/src/language/modules/DevUtilsModule.cpp
+++ b/src/language/modules/DevUtilsModule.cpp
@@ -24,6 +24,16 @@ template <>
 inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemArrayVariant>> =
   ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_array");
 
+class SubItemValuePerItemVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const SubItemValuePerItemVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("sub_item_value");
+
+class SubItemArrayPerItemVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const SubItemArrayPerItemVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("sub_item_array");
+
 DevUtilsModule::DevUtilsModule()
 {
   this->_addBuiltinFunction("getAST", std::function(
@@ -97,6 +107,37 @@ DevUtilsModule::DevUtilsModule()
                                 parallel_check(*item_value, name, ASTBacktrace::getInstance().sourceLocation());
                               }
 
+                              ));
+
+  this->_addBuiltinFunction("parallel_check",
+                            std::function(
+
+                              [](const std::shared_ptr<const ItemArrayVariant>& item_array, const std::string& name) {
+                                parallel_check(*item_array, name, ASTBacktrace::getInstance().sourceLocation());
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("parallel_check",
+                            std::function(
+
+                              [](const std::shared_ptr<const SubItemValuePerItemVariant>& subitem_value_per_item,
+                                 const std::string& name) {
+                                parallel_check(*subitem_value_per_item, name,
+                                               ASTBacktrace::getInstance().sourceLocation());
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("parallel_check",
+                            std::function(
+
+                              [](const std::shared_ptr<const SubItemArrayPerItemVariant>& subitem_array_per_item,
+                                 const std::string& name) {
+                                parallel_check(*subitem_array_per_item, name,
+                                               ASTBacktrace::getInstance().sourceLocation());
+                              }
+
                               ));
 }
 
diff --git a/src/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index a98c039e0..14e0caa63 100644
--- a/src/mesh/SubItemArrayPerItem.hpp
+++ b/src/mesh/SubItemArrayPerItem.hpp
@@ -49,6 +49,13 @@ class SubItemArrayPerItem
   friend SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
  public:
+  // This is not the correct way to look at SubItemArrayPerItem, use with care
+  Table<const DataType>
+  tableView() const
+  {
+    return m_values;
+  }
+
   [[nodiscard]] friend PUGS_INLINE SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
diff --git a/src/mesh/SubItemArrayPerItemVariant.hpp b/src/mesh/SubItemArrayPerItemVariant.hpp
index ffb0ea673..2f917aa7c 100644
--- a/src/mesh/SubItemArrayPerItemVariant.hpp
+++ b/src/mesh/SubItemArrayPerItemVariant.hpp
@@ -146,7 +146,7 @@ class SubItemArrayPerItemVariant
  public:
   PUGS_INLINE
   const Variant&
-  itemArray() const
+  subItemArrayPerItem() const
   {
     return m_sub_item_array_per_item;
   }
@@ -189,7 +189,7 @@ class SubItemArrayPerItemVariant
                   "SubItemArrayPerItem with this DataType is not allowed in variant");
   }
 
-  SubItemArrayPerItemVariant& operator=(SubItemArrayPerItemVariant&&) = default;
+  SubItemArrayPerItemVariant& operator=(SubItemArrayPerItemVariant&&)      = default;
   SubItemArrayPerItemVariant& operator=(const SubItemArrayPerItemVariant&) = default;
 
   SubItemArrayPerItemVariant(const SubItemArrayPerItemVariant&) = default;
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index f5ce39d5f..d1b569919 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -48,6 +48,14 @@ class SubItemValuePerItem
   friend SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
  public:
+  // This is not the correct way to look at SubItemValuePerItem, use
+  // with care
+  Array<const DataType>
+  arrayView() const
+  {
+    return m_values;
+  }
+
   [[nodiscard]] friend PUGS_INLINE SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
diff --git a/src/mesh/SubItemValuePerItemVariant.hpp b/src/mesh/SubItemValuePerItemVariant.hpp
index 618fd82d9..ef2bf58e3 100644
--- a/src/mesh/SubItemValuePerItemVariant.hpp
+++ b/src/mesh/SubItemValuePerItemVariant.hpp
@@ -146,7 +146,7 @@ class SubItemValuePerItemVariant
  public:
   PUGS_INLINE
   const Variant&
-  itemValue() const
+  subItemValuePerItem() const
   {
     return m_sub_item_value_per_item;
   }
@@ -189,7 +189,7 @@ class SubItemValuePerItemVariant
                   "SubItemValuePerItem with this DataType is not allowed in variant");
   }
 
-  SubItemValuePerItemVariant& operator=(SubItemValuePerItemVariant&&) = default;
+  SubItemValuePerItemVariant& operator=(SubItemValuePerItemVariant&&)      = default;
   SubItemValuePerItemVariant& operator=(const SubItemValuePerItemVariant&) = default;
 
   SubItemValuePerItemVariant(const SubItemValuePerItemVariant&) = default;
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index 54bec14a0..104e70cdd 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -120,8 +120,8 @@ initialize(int& argc, char* argv[])
     app.add_flag("--exec-stat,!--no-exec-stat", print_exec_stat,
                  "Display memory and CPU usage after execution [default: true]");
 
-    bool show_backtrace = true;
-    app.add_flag("-b,--backtrace,!--no-backtrace", show_backtrace, "Show backtrace on failure [default: true]");
+    bool show_backtrace = false;
+    app.add_flag("-b,--backtrace,!--no-backtrace", show_backtrace, "Show backtrace on failure [default: false]");
 
     app.add_flag("--signal,!--no-signal", enable_signals, "Catches signals [default: true]");
 
-- 
GitLab