diff --git a/CMakeLists.txt b/CMakeLists.txt
index ee35378147881930584a2c58444c9a1d4f4d64ff..b7f6f72f6cbef2278363d53cfd1f02ee04758421 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,7 +14,7 @@ find_package(PkgConfig REQUIRED)
 #------------------------------------------------------
 
 # custom variable allowing to define version suffixes such as -rc*, -beta*, ...
-set(PUGS_VERSION "0.4.1")
+set(PUGS_VERSION "0.5.0")
 
 # deduce PUGS_SHORT_VERSION using regex
 string(REGEX MATCH "^[0-9]+\.[0-9]+\.[0-9]+" PUGS_SHORT_VERSION ${PUGS_VERSION})
diff --git a/src/language/PugsParser.cpp b/src/language/PugsParser.cpp
index cb42fc45acf4f4bfaf14fe2c11a3b0ab84938af6..fb7f1926cae6c93cb9d0bbfed7435271133d6a6d 100644
--- a/src/language/PugsParser.cpp
+++ b/src/language/PugsParser.cpp
@@ -42,6 +42,15 @@
 #include <unordered_map>
 #include <variant>
 
+void
+clear_symbol_tables(std::unique_ptr<ASTNode>& node)
+{
+  for (auto& child : node->children) {
+    clear_symbol_tables(child);
+  }
+  node->m_symbol_table->clearValues();
+}
+
 void
 parser(const std::string& filename)
 {
@@ -103,6 +112,7 @@ parser(const std::string& filename)
     }
     catch (language::Exit& e) {
       ExecutionStatManager::getInstance().setExitCode(e.code());
+      clear_symbol_tables(root_node);
     }
 
     root_node->m_symbol_table->clearValues();
diff --git a/src/mesh/ItemArrayUtils.hpp b/src/mesh/ItemArrayUtils.hpp
index 0d90f77af27add5f08982e2be535bd7b779837a9..cb8ee339b56a3f127b3f446b7ee28970b31aa143 100644
--- a/src/mesh/ItemArrayUtils.hpp
+++ b/src/mesh/ItemArrayUtils.hpp
@@ -165,7 +165,7 @@ max(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array)
     void
     init(data_type& value) const
     {
-      value = std::numeric_limits<data_type>::min();
+      value = std::numeric_limits<data_type>::lowest();
     }
 
     PUGS_INLINE
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index e64f2e1bdc0bc1e934aac22258461947c1c0ffd3..031caa640dbc8617ca4514c513035975451f7939 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -156,7 +156,7 @@ max(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
     void
     init(data_type& value) const
     {
-      value = std::numeric_limits<data_type>::min();
+      value = std::numeric_limits<data_type>::lowest();
     }
 
     PUGS_INLINE
diff --git a/src/mesh/MeshFlatEdgeBoundary.cpp b/src/mesh/MeshFlatEdgeBoundary.cpp
index c20b9c1a5967fa511fdd4c9f99359696d7117cc2..7a2d4918d9c822303ba6debee91b5f4d707569c4 100644
--- a/src/mesh/MeshFlatEdgeBoundary.cpp
+++ b/src/mesh/MeshFlatEdgeBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshFlatEdgeBoundary<MeshType>
 getMeshFlatEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatEdgeBoundary.hpp b/src/mesh/MeshFlatEdgeBoundary.hpp
index 10ab8edb835e8c7652665cae4b9db9cfe75c4418..70b51bc982222b7006b21e3a825cfeeab3fe6c93 100644
--- a/src/mesh/MeshFlatEdgeBoundary.hpp
+++ b/src/mesh/MeshFlatEdgeBoundary.hpp
@@ -24,7 +24,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   MeshFlatEdgeBoundary& operator=(const MeshFlatEdgeBoundary&) = default;
   MeshFlatEdgeBoundary& operator=(MeshFlatEdgeBoundary&&)      = default;
 
-  template <MeshConcept MeshTypeT>
+  template <typename MeshTypeT>
   friend MeshFlatEdgeBoundary<MeshTypeT> getMeshFlatEdgeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -40,7 +40,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   ~MeshFlatEdgeBoundary()                           = default;
 };
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshFlatEdgeBoundary<MeshType> getMeshFlatEdgeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshFlatFaceBoundary.cpp b/src/mesh/MeshFlatFaceBoundary.cpp
index 632c6a033dee7fc5593363660cc179faa966fb80..2ecdf6768eb865413386b10a1a419c6e30799def 100644
--- a/src/mesh/MeshFlatFaceBoundary.cpp
+++ b/src/mesh/MeshFlatFaceBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshFlatFaceBoundary<MeshType>
 getMeshFlatFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatFaceBoundary.hpp b/src/mesh/MeshFlatFaceBoundary.hpp
index 0c53a31e095251c6ffaa70b7f6d58bfce10a6864..fc4d9d0f82b15a0d02dea4efdff3a13a6cceeb7a 100644
--- a/src/mesh/MeshFlatFaceBoundary.hpp
+++ b/src/mesh/MeshFlatFaceBoundary.hpp
@@ -24,7 +24,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   MeshFlatFaceBoundary& operator=(const MeshFlatFaceBoundary&) = default;
   MeshFlatFaceBoundary& operator=(MeshFlatFaceBoundary&&)      = default;
 
-  template <MeshConcept MeshTypeT>
+  template <typename MeshTypeT>
   friend MeshFlatFaceBoundary<MeshTypeT> getMeshFlatFaceBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -40,7 +40,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   ~MeshFlatFaceBoundary()                           = default;
 };
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshFlatFaceBoundary<MeshType> getMeshFlatFaceBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index bdc7c20a2c854c2a9a6492d8f06f46a026b38b0b..c3ca1dbd4ce980e6290b99af8a7bca878ca4520f 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -323,7 +323,7 @@ MeshFlatNodeBoundary<Mesh<3>>::_getOutgoingNormal(const Mesh<3>& mesh)
   }
 }
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshFlatNodeBoundary<MeshType>
 getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
index 2eb28d14899a1d70f47f8cd139b435c4a81390b7..66c610ea72a3aed5e6876516890c0d168bc9aac4 100644
--- a/src/mesh/MeshFlatNodeBoundary.hpp
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -20,21 +20,23 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
   Rd _getNormal(const MeshType& mesh);
 
   void _checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
-                            const TinyVector<MeshType::Dimension, double>& origin, const double length,
+                            const TinyVector<MeshType::Dimension, double>& origin,
+                            const double length,
                             const MeshType& mesh) const;
 
   Rd _getOutgoingNormal(const MeshType& mesh);
 
  public:
-  const Rd& outgoingNormal() const
+  const Rd&
+  outgoingNormal() const
   {
     return m_outgoing_normal;
   }
 
   MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
+  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&)      = default;
 
-  template <MeshConcept MeshTypeT>
+  template <typename MeshTypeT>
   friend MeshFlatNodeBoundary<MeshTypeT> getMeshFlatNodeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -50,11 +52,11 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
  public:
   MeshFlatNodeBoundary()                            = default;
   MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary(MeshFlatNodeBoundary &&)     = default;
+  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
   ~MeshFlatNodeBoundary()                           = default;
 };
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshFlatNodeBoundary<MeshType> getMeshFlatNodeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineEdgeBoundary.cpp b/src/mesh/MeshLineEdgeBoundary.cpp
index 40de23d1e2dc86bbd82b8542b75e32d1cac40352..f1e0333a9afd8a642f3809560f0c16431e98e1e4 100644
--- a/src/mesh/MeshLineEdgeBoundary.cpp
+++ b/src/mesh/MeshLineEdgeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshLineEdgeBoundary<MeshType>
 getMeshLineEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineEdgeBoundary.hpp b/src/mesh/MeshLineEdgeBoundary.hpp
index 0fc5db0ba6b73960cd1a64c5a0626585c857a4b8..11994eefa1d84580d6c21b58a4d5332ba8479431 100644
--- a/src/mesh/MeshLineEdgeBoundary.hpp
+++ b/src/mesh/MeshLineEdgeBoundary.hpp
@@ -18,7 +18,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
   const Rd m_direction;
 
  public:
-  template <MeshConcept MeshTypeT>
+  template <typename MeshTypeT>
   friend MeshLineEdgeBoundary<MeshTypeT> getMeshLineEdgeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -44,7 +44,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
   ~MeshLineEdgeBoundary()                           = default;
 };
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshLineEdgeBoundary<MeshType> getMeshLineEdgeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineFaceBoundary.cpp b/src/mesh/MeshLineFaceBoundary.cpp
index 3608e7aef37dfeddf9eb2c6ad7d207d1acae3f07..81ea991ed6dcf999e3144a9b4d9d48f080a1499e 100644
--- a/src/mesh/MeshLineFaceBoundary.cpp
+++ b/src/mesh/MeshLineFaceBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshLineFaceBoundary<MeshType>
 getMeshLineFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineFaceBoundary.hpp b/src/mesh/MeshLineFaceBoundary.hpp
index 206961a08f88f9a4008a14a7a223b9d90ff323ed..5b993923c2836500f46ebbc112ef234c05c67f94 100644
--- a/src/mesh/MeshLineFaceBoundary.hpp
+++ b/src/mesh/MeshLineFaceBoundary.hpp
@@ -17,7 +17,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
   const Rd m_direction;
 
  public:
-  template <MeshConcept MeshTypeT>
+  template <typename MeshTypeT>
   friend MeshLineFaceBoundary<MeshTypeT> getMeshLineFaceBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -43,7 +43,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
   ~MeshLineFaceBoundary()                           = default;
 };
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshLineFaceBoundary<MeshType> getMeshLineFaceBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
index 1cd3bb64d749d49de82d925bff9eaa3ecf88ead9..7ec8bb898aa7af4a8d747e516f0438b1792edb67 100644
--- a/src/mesh/MeshLineNodeBoundary.cpp
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -105,7 +105,7 @@ MeshLineNodeBoundary<Mesh<3>>::_getDirection(const Mesh<3>& mesh)
   return direction;
 }
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshLineNodeBoundary<MeshType>
 getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
index fb75ef9a49d77095c47d1a45db42e636b905c572..95cbab34cfbb35190fe23f607bdbcd46698f4023 100644
--- a/src/mesh/MeshLineNodeBoundary.hpp
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -24,7 +24,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
                             const MeshType& mesh) const;
 
  public:
-  template <MeshConcept MeshTypeT>
+  template <typename MeshTypeT>
   friend MeshLineNodeBoundary<MeshTypeT> getMeshLineNodeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -58,7 +58,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
   ~MeshLineNodeBoundary()                           = default;
 };
 
-template <MeshConcept MeshType>
+template <typename MeshType>
 MeshLineNodeBoundary<MeshType> getMeshLineNodeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshNodeBoundaryUtils.cpp b/src/mesh/MeshNodeBoundaryUtils.cpp
index d6a0e31cd9a671530216b9dce15747e07f06c9b4..39715d33d2f80da11dc9b11d1433320305fe2715 100644
--- a/src/mesh/MeshNodeBoundaryUtils.cpp
+++ b/src/mesh/MeshNodeBoundaryUtils.cpp
@@ -17,7 +17,7 @@ getBounds(const Mesh<2>& mesh, const RefNodeList& ref_node_list)
   R2& xmax = bounds[1];
 
   xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
+  xmax = R2{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest()};
 
   auto update_xmin = [](const R2& x, R2& x_min) {
     if ((x[0] < x_min[0]) or ((x[0] == x_min[0]) and (x[1] < x_min[1]))) {
@@ -118,12 +118,12 @@ getBounds(const Mesh<3>& mesh, const RefNodeList& ref_node_list)
   ymin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
   zmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
 
-  xmax =
-    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  ymax =
-    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  zmax =
-    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  xmax = R3{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
+            std::numeric_limits<double>::lowest()};
+  ymax = R3{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
+            std::numeric_limits<double>::lowest()};
+  zmax = R3{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
+            std::numeric_limits<double>::lowest()};
 
   const NodeValue<const R3>& xr = mesh.xr();
 
diff --git a/src/mesh/SubItemArrayPerItemUtils.hpp b/src/mesh/SubItemArrayPerItemUtils.hpp
index aec08ef845cc80bafd20ec8fa1a3b9bcc94415de..70b835ae3d8ff75c0fad7c9cd1b13c6cfb8ad998 100644
--- a/src/mesh/SubItemArrayPerItemUtils.hpp
+++ b/src/mesh/SubItemArrayPerItemUtils.hpp
@@ -170,7 +170,7 @@ max(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_a
     void
     init(data_type& value) const
     {
-      value = std::numeric_limits<data_type>::min();
+      value = std::numeric_limits<data_type>::lowest();
     }
 
     PUGS_INLINE
diff --git a/src/mesh/SubItemValuePerItemUtils.hpp b/src/mesh/SubItemValuePerItemUtils.hpp
index 575c39c5238877bb064883c54b7d08792547000b..cc0c941f7f23ee5592e861d2a82a23ca07b3d0e1 100644
--- a/src/mesh/SubItemValuePerItemUtils.hpp
+++ b/src/mesh/SubItemValuePerItemUtils.hpp
@@ -166,7 +166,7 @@ max(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_v
     void
     init(data_type& value) const
     {
-      value = std::numeric_limits<data_type>::min();
+      value = std::numeric_limits<data_type>::lowest();
     }
 
     PUGS_INLINE
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index cd7834d713d2f9ca94d6dd65aa301ed444c6340f..d33c5e0d2b6f2eaeb5af453ce3c388336f03c61c 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -222,7 +222,7 @@ struct VTKWriter::VTKType
       if constexpr (std::is_unsigned_v<DataType>) {
         return "UInt" + stringify(sizeof(DataType) * 8);
       } else {
-        return "UInt" + stringify(sizeof(DataType) * 8);
+        return "Int" + stringify(sizeof(DataType) * 8);
       }
     } else if constexpr (std::is_floating_point_v<DataType>) {
       return "Float" + stringify(sizeof(DataType) * 8);
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
index 1964068a0222bfb35304d30afab2861f3cb97383..c07c20c711e827ae566635ca82c9456bbda853ac 100644
--- a/src/output/WriterBase.hpp
+++ b/src/output/WriterBase.hpp
@@ -61,7 +61,7 @@ class WriterBase : public IWriter
       if (m_saved_times.size() > 0) {
         return m_saved_times[m_saved_times.size() - 1];
       } else {
-        return -std::numeric_limits<double>::max();
+        return std::numeric_limits<double>::lowest();
       }
     }
 
@@ -77,7 +77,8 @@ class WriterBase : public IWriter
 
     PeriodManager(const PeriodManager&) = default;
     PeriodManager(PeriodManager&&)      = default;
-    PeriodManager(double time_period) : m_time_period{time_period}, m_next_time{-std::numeric_limits<double>::max()} {}
+    PeriodManager(double time_period) : m_time_period{time_period}, m_next_time{std::numeric_limits<double>::lowest()}
+    {}
   };
 
  protected:
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 90998d3aba2e9ecf0a9b5c4e8961c7b87889eaa5..e4db63bd7ec57502c08a1e3f212386ebe8215859 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -371,8 +371,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                  const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
-    if (not i_mesh) {
+    std::shared_ptr mesh_v = getCommonMesh({rho_v, c_v, u_v, p_v});
+    if (not mesh_v) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
@@ -380,7 +380,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
       throw NormalError("acoustic solver expects P0 functions");
     }
 
-    const MeshType& mesh              = *i_mesh->get<MeshType>();
+    const MeshType& mesh              = *mesh_v->get<MeshType>();
     const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
     const DiscreteScalarFunction& c   = c_v->get<DiscreteScalarFunction>();
     const DiscreteVectorFunction& u   = u_v->get<DiscreteVectorFunction>();
@@ -474,8 +474,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
-    if (not i_mesh) {
+    std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v});
+    if (not mesh_v) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
@@ -484,7 +484,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     }
 
     return this->apply_fluxes(dt,                                     //
-                              *i_mesh->get<MeshType>(),               //
+                              *mesh_v->get<MeshType>(),               //
                               rho_v->get<DiscreteScalarFunction>(),   //
                               u_v->get<DiscreteVectorFunction>(),     //
                               E_v->get<DiscreteScalarFunction>(),     //
@@ -863,9 +863,9 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::SymmetryBoundaryCondition
   ~SymmetryBoundaryCondition() = default;
 };
 
-AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const MeshVariant>& i_mesh)
+AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v)
 {
-  if (not i_mesh) {
+  if (not mesh_v) {
     throw NormalError("discrete functions are not defined on the same mesh");
   }
 
@@ -878,5 +878,5 @@ AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const MeshVar
         throw NormalError("unexpected mesh type");
       }
     },
-    i_mesh->variant());
+    mesh_v->variant());
 }
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index d88067f982e244a53cb6f87d062e404392188a69..7bc6cd8c0dc854ecaa8caa2d8809891ed328790d 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -388,8 +388,8 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
-    if (not i_mesh) {
+    std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+    if (not mesh_v) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
@@ -397,7 +397,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
       throw NormalError("hyperelastic solver expects P0 functions");
     }
 
-    const MeshType& mesh                = *i_mesh->get<MeshType>();
+    const MeshType& mesh                = *mesh_v->get<MeshType>();
     const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
     const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
     const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
@@ -505,8 +505,8 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
-    if (not i_mesh) {
+    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
+    if (not mesh_v) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
@@ -515,7 +515,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
     }
 
     return this->apply_fluxes(dt,                                   //
-                              *i_mesh->get<MeshType>(),             //
+                              *mesh_v->get<MeshType>(),             //
                               rho->get<DiscreteScalarFunction>(),   //
                               u->get<DiscreteVectorFunction>(),     //
                               E->get<DiscreteScalarFunction>(),     //
@@ -939,9 +939,9 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::SymmetryBoundaryC
   ~SymmetryBoundaryCondition() = default;
 };
 
-HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& i_mesh)
+HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v)
 {
-  if (not i_mesh) {
+  if (not mesh_v) {
     throw NormalError("discrete functions are not defined on the same mesh");
   }
 
@@ -954,5 +954,5 @@ HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const
         throw NormalError("unexpected mesh type");
       }
     },
-    i_mesh->variant());
+    mesh_v->variant());
 }
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 2f0146e2b9d08c231eaadee7e7f50eeea439ddaf..a08453ba766450909f3bcdf0f0fd678bdb06456c 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -355,7 +355,7 @@ max(const Array<DataType>& array)
     void
     init(data_type& value) const
     {
-      value = std::numeric_limits<data_type>::min();
+      value = std::numeric_limits<data_type>::lowest();
     }
 
     PUGS_INLINE
diff --git a/src/utils/SmallArray.hpp b/src/utils/SmallArray.hpp
index 9e0782f1c0adb0b9bfa4bcf25bcaa8977a3800d4..81d70a848187216e82cf80c5208f96cb64e05dfa 100644
--- a/src/utils/SmallArray.hpp
+++ b/src/utils/SmallArray.hpp
@@ -25,12 +25,14 @@ class [[nodiscard]] SmallArray
   friend SmallArray<std::add_const_t<DataType>>;
 
  public:
-  PUGS_INLINE size_t size() const noexcept
+  PUGS_INLINE size_t
+  size() const noexcept
   {
     return m_size;
   }
 
-  friend PUGS_INLINE SmallArray<std::remove_const_t<DataType>> copy(const SmallArray<DataType>& source)
+  friend PUGS_INLINE SmallArray<std::remove_const_t<DataType>>
+  copy(const SmallArray<DataType>& source)
   {
     SmallArray<std::remove_const_t<DataType>> image(source.m_size);
     std::copy(source.m_values.get(), source.m_values.get() + source.m_size, image.m_values.get());
@@ -38,28 +40,31 @@ class [[nodiscard]] SmallArray
     return image;
   }
 
-  friend PUGS_INLINE void copy_to(const SmallArray<DataType>& source,
-                                  const SmallArray<std::remove_const_t<DataType>>& destination)
+  friend PUGS_INLINE void
+  copy_to(const SmallArray<DataType>& source, const SmallArray<std::remove_const_t<DataType>>& destination)
   {
     Assert(source.size() == destination.size());
     std::copy(source.m_values.get(), source.m_values.get() + source.m_size, destination.m_values.get());
   }
 
-  PUGS_INLINE DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  PUGS_INLINE DataType&
+  operator[](index_type i) const noexcept(NO_ASSERT)
   {
     Assert(i < m_size);
     return m_values[i];
   }
 
   PUGS_INLINE
-  void fill(const DataType& data) const
+  void
+  fill(const DataType& data) const
   {
     static_assert(not std::is_const_v<DataType>, "Cannot modify SmallArray of const");
     std::fill(m_values.get(), m_values.get() + m_size, data);
   }
 
   template <typename DataType2>
-  PUGS_INLINE SmallArray& operator=(const SmallArray<DataType2>& array) noexcept
+  PUGS_INLINE SmallArray&
+  operator=(const SmallArray<DataType2>& array) noexcept
   {
     // ensures that DataType is the same as source DataType2
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
@@ -99,7 +104,8 @@ class [[nodiscard]] SmallArray
 #endif   // NDEBUG
   }
 
-  friend std::ostream& operator<<(std::ostream& os, const SmallArray& x)
+  friend std::ostream&
+  operator<<(std::ostream& os, const SmallArray& x)
   {
     if (x.size() > 0) {
       os << 0 << ':' << NaNHelper(x[0]);
@@ -117,13 +123,14 @@ class [[nodiscard]] SmallArray
   SmallArray(const SmallArray&) = default;
 
   template <typename DataType2>
-  PUGS_INLINE SmallArray(const SmallArray<DataType2>& array) noexcept
+  PUGS_INLINE
+  SmallArray(const SmallArray<DataType2>& array) noexcept
   {
     this->operator=(array);
   }
 
   PUGS_INLINE
-  SmallArray(SmallArray &&) = default;
+  SmallArray(SmallArray&&) = default;
 
   PUGS_INLINE
   ~SmallArray() = default;
@@ -172,7 +179,7 @@ max(const SmallArray<DataType>& array)
   using data_type  = std::remove_const_t<DataType>;
   using index_type = typename SmallArray<DataType>::index_type;
 
-  data_type max_value = -std::numeric_limits<data_type>::max();
+  data_type max_value = std::numeric_limits<data_type>::lowest();
   for (index_type i = 0; i < array.size(); ++i) {
     if (array[i] > max_value) {
       max_value = array[i];
diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
index 4a447f70a9e99151027266602763b712d00525cf..28c018cced7ac2c4a2295b1824882a7dc5b492f8 100644
--- a/tests/test_Array.cpp
+++ b/tests/test_Array.cpp
@@ -324,6 +324,103 @@ TEST_CASE("Array", "[utils]")
     }
   }
 
+  SECTION("checking for floating Array min/max")
+  {
+    SECTION("Min")
+    {
+      Array<double> b(10);
+      b[0] = 13;
+      b[1] = 1;
+      b[2] = 8;
+      b[3] = -3;
+      b[4] = 23;
+      b[5] = -1;
+      b[6] = 13;
+      b[7] = 0;
+      b[8] = 12;
+      b[9] = 9;
+
+      REQUIRE(min(b) == -3);
+
+      b.fill(0);
+      REQUIRE(min(b) == 0);
+
+      b[0] = -13;
+      b[1] = -1;
+      b[2] = -8;
+      b[3] = -3;
+      b[4] = -23;
+      b[5] = -1;
+      b[6] = -13;
+      b[7] = 0;
+      b[8] = -12;
+      b[9] = -9;
+
+      REQUIRE(min(b) == -23);
+
+      b[0] = 13;
+      b[1] = 1;
+      b[2] = 8;
+      b[3] = 3;
+      b[4] = 23;
+      b[5] = 1;
+      b[6] = 13;
+      b[7] = 0;
+      b[8] = 12;
+      b[9] = 9;
+
+      REQUIRE(min(b) == 0);
+
+      b[7] = 3;
+      REQUIRE(min(b) == 1);
+    }
+
+    SECTION("Max")
+    {
+      Array<double> b(10);
+      b[0] = 13;
+      b[1] = 1;
+      b[2] = 8;
+      b[3] = -3;
+      b[4] = 23;
+      b[5] = -1;
+      b[6] = 13;
+      b[7] = 0;
+      b[8] = 12;
+      b[9] = 9;
+
+      REQUIRE(max(b) == 23);
+
+      b[0] = -13;
+      b[1] = -12;
+      b[2] = -8;
+      b[3] = -3;
+      b[4] = -23;
+      b[5] = -1;
+      b[6] = -13;
+      b[7] = -10;
+      b[8] = -12;
+      b[9] = -9;
+
+      REQUIRE(max(b) == -1);
+
+      b.fill(-13);
+      REQUIRE(max(b) == -13);
+
+      b[0] = 13;
+      b[1] = 12;
+      b[2] = 8;
+      b[3] = 3;
+      b[4] = 23;
+      b[5] = 1;
+      b[6] = 13;
+      b[7] = 10;
+      b[8] = 12;
+      b[9] = 9;
+      REQUIRE(max(b) == 23);
+    }
+  }
+
   SECTION("reproducible floating point sum")
   {
     auto direct_sum = [](auto array) {
diff --git a/tests/test_ItemArrayUtils.cpp b/tests/test_ItemArrayUtils.cpp
index 8d04842849ed469e1cf0c1e6e217d8e1da4f01f8..0237b90da16549c30fd8b581643caec3dc287ddd 100644
--- a/tests/test_ItemArrayUtils.cpp
+++ b/tests/test_ItemArrayUtils.cpp
@@ -1028,6 +1028,66 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         }
       }
     }
+
+    SECTION("max for negative double values")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d_v = named_mesh.mesh();
+          auto mesh_3d   = mesh_3d_v->get<Mesh<3>>();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          CellArray<double> cell_array{connectivity, 2};
+          cell_array.fill(std::numeric_limits<double>::max());
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            mesh_3d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+                  cell_array[cell_id][i] = -parallel::rank() - 10 + i;
+                }
+              }
+            });
+
+          REQUIRE(max(cell_array) == -parallel::size() - 10);
+        }
+      }
+    }
+
+    SECTION("max for negative integral values")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d_v = named_mesh.mesh();
+          auto mesh_3d   = mesh_3d_v->get<Mesh<3>>();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          CellArray<int> cell_array{connectivity, 2};
+          cell_array.fill(std::numeric_limits<int>::max());
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            mesh_3d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+                  cell_array[cell_id][i] = -2 * parallel::size() + parallel::rank() - 10 + i;
+                }
+              }
+            });
+
+          REQUIRE(max(cell_array) == -static_cast<int>(parallel::size() + 10));
+        }
+      }
+    }
   }
 
   SECTION("sum")
diff --git a/tests/test_ItemValueUtils.cpp b/tests/test_ItemValueUtils.cpp
index f80819f2c619709584c6574dd185d5a0e2f512fb..149dcbfa6c7c6d7976b98f21b559a03016daf857 100644
--- a/tests/test_ItemValueUtils.cpp
+++ b/tests/test_ItemValueUtils.cpp
@@ -240,6 +240,34 @@ TEST_CASE("ItemValueUtils", "[mesh]")
         }
       }
     }
+
+    SECTION("max of all negative values")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d_v = named_mesh.mesh();
+          auto mesh_3d   = mesh_3d_v->get<Mesh<3>>();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          CellValue<float> cell_value{connectivity};
+          cell_value.fill(std::numeric_limits<float>::max());
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            mesh_3d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                cell_value[cell_id] = -2 * parallel::size() + parallel::rank() + 1;
+              }
+            });
+
+          REQUIRE(max(cell_value) == -parallel::size());
+        }
+      }
+    }
   }
 
   SECTION("sum")
diff --git a/tests/test_SubItemArrayPerItemUtils.cpp b/tests/test_SubItemArrayPerItemUtils.cpp
index 1ff8bb17f7bf170c748cfa15aa3fc82f25f4a879..ce3bb6265fd3e48a13d0505e39ddfa312bf33e2f 100644
--- a/tests/test_SubItemArrayPerItemUtils.cpp
+++ b/tests/test_SubItemArrayPerItemUtils.cpp
@@ -283,6 +283,39 @@ TEST_CASE("SubItemArrayPerItemUtils", "[mesh]")
         }
       }
     }
+
+    SECTION("max of only negative values")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh()->get<Mesh<3>>();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          NodeArrayPerEdge<double> node_array_per_edge{connectivity, 3};
+          node_array_per_edge.fill(std::numeric_limits<double>::max());
+
+          auto edge_is_owned = connectivity.edgeIsOwned();
+
+          parallel_for(
+            connectivity.numberOfEdges(), PUGS_LAMBDA(EdgeId edge_id) {
+              if (edge_is_owned[edge_id]) {
+                auto edge_table = node_array_per_edge.itemTable(edge_id);
+                for (size_t i = 0; i < edge_table.numberOfRows(); ++i) {
+                  for (size_t j = 0; j < edge_table.numberOfColumns(); ++j) {
+                    edge_table(i, j) = -10. - parallel::size() + parallel::rank() - i - j;
+                  }
+                }
+              }
+            });
+
+          REQUIRE(max(node_array_per_edge) == -11.);
+        }
+      }
+    }
   }
 
   SECTION("sum")
diff --git a/tests/test_SubItemValuePerItemUtils.cpp b/tests/test_SubItemValuePerItemUtils.cpp
index 69214c1c933ba6abb0c93e094317901263e70eed..fd84ad62e9dcfb1f5dd379c8f469a5b5d7e04d0d 100644
--- a/tests/test_SubItemValuePerItemUtils.cpp
+++ b/tests/test_SubItemValuePerItemUtils.cpp
@@ -267,6 +267,36 @@ TEST_CASE("SubItemValuePerItemUtils", "[mesh]")
         }
       }
     }
+
+    SECTION("max for all negative values")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh()->get<Mesh<2>>();
+
+          const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+          EdgeValuePerCell<double> edge_value_per_cell{connectivity};
+          edge_value_per_cell.fill(std::numeric_limits<double>::max());
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                auto cell_array = edge_value_per_cell.itemArray(cell_id);
+                for (size_t i = 0; i < cell_array.size(); ++i) {
+                  cell_array[i] = -1. * parallel::size() - 10 + parallel::rank() - i;
+                }
+              }
+            });
+
+          REQUIRE(max(edge_value_per_cell) == -11);
+        }
+      }
+    }
   }
 
   SECTION("sum")