diff --git a/CMakeLists.txt b/CMakeLists.txt
index ef624927afc0f23aeec959727bb1bec3ed627674..d7c7318458823c05c9f385c1b5c50118534f6a76 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -257,8 +257,8 @@ else ()
 endif()
 
 #------------------------------------------------------
-# C++ 17 flags
-set(CMAKE_CXX_STANDARD "17")
+# C++ 20 flags
+set(CMAKE_CXX_STANDARD "20")
 
 #------------------------------------------------------
 # Kokkos configuration
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 8e0bbd4db39f20977afa3725177f728d7094710d..94590ff656784d018e9ca57991674da773d64568 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -262,15 +262,6 @@ class [[nodiscard]] TinyMatrix
     return *this;
   }
 
-  PUGS_INLINE
-  constexpr void
-  operator+=(const volatile TinyMatrix& A) volatile
-  {
-    for (size_t i = 0; i < M * N; ++i) {
-      m_values[i] += A.m_values[i];
-    }
-  }
-
   PUGS_INLINE
   constexpr TinyMatrix&
   operator-=(const TinyMatrix& A)
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index e39dcec6e9834e5cb5eb80e1e845de1054beb1b2..2bc9c79c36ebe2caef6cf63ac6dc653d3a8b0bcb 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -23,7 +23,8 @@ class [[nodiscard]] TinyVector
   static_assert((N > 0), "TinyVector size must be strictly positive");
 
   template <typename... Args>
-  PUGS_FORCEINLINE constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
+  PUGS_FORCEINLINE constexpr void
+  _unpackVariadicInput(const T& t, Args&&... args) noexcept
   {
     m_values[N - 1 - sizeof...(args)] = t;
     if constexpr (sizeof...(args) > 0) {
@@ -33,7 +34,8 @@ class [[nodiscard]] TinyVector
 
  public:
   PUGS_INLINE
-  constexpr TinyVector operator-() const
+  constexpr TinyVector
+  operator-() const
   {
     TinyVector opposite;
     for (size_t i = 0; i < N; ++i) {
@@ -42,12 +44,14 @@ class [[nodiscard]] TinyVector
     return opposite;
   }
 
-  [[nodiscard]] PUGS_INLINE constexpr size_t dimension() const
+  [[nodiscard]] PUGS_INLINE constexpr size_t
+  dimension() const
   {
     return N;
   }
 
-  [[nodiscard]] PUGS_INLINE constexpr bool operator==(const TinyVector& v) const
+  [[nodiscard]] PUGS_INLINE constexpr bool
+  operator==(const TinyVector& v) const
   {
     for (size_t i = 0; i < N; ++i) {
       if (m_values[i] != v.m_values[i])
@@ -56,12 +60,14 @@ class [[nodiscard]] TinyVector
     return true;
   }
 
-  [[nodiscard]] PUGS_INLINE constexpr bool operator!=(const TinyVector& v) const
+  [[nodiscard]] PUGS_INLINE constexpr bool
+  operator!=(const TinyVector& v) const
   {
     return not this->operator==(v);
   }
 
-  [[nodiscard]] PUGS_INLINE constexpr friend T dot(const TinyVector& u, const TinyVector& v)
+  [[nodiscard]] PUGS_INLINE constexpr friend T
+  dot(const TinyVector& u, const TinyVector& v)
   {
     T t = u.m_values[0] * v.m_values[0];
     for (size_t i = 1; i < N; ++i) {
@@ -71,7 +77,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr TinyVector& operator*=(const T& t)
+  constexpr TinyVector&
+  operator*=(const T& t)
   {
     for (size_t i = 0; i < N; ++i) {
       m_values[i] *= t;
@@ -80,21 +87,24 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr friend TinyVector operator*(const T& t, const TinyVector& v)
+  constexpr friend TinyVector
+  operator*(const T& t, const TinyVector& v)
   {
     TinyVector w = v;
     return w *= t;
   }
 
   PUGS_INLINE
-  constexpr friend TinyVector operator*(const T& t, TinyVector&& v)
+  constexpr friend TinyVector
+  operator*(const T& t, TinyVector&& v)
   {
     v *= t;
     return std::move(v);
   }
 
   PUGS_INLINE
-  constexpr friend std::ostream& operator<<(std::ostream& os, const TinyVector& v)
+  constexpr friend std::ostream&
+  operator<<(std::ostream& os, const TinyVector& v)
   {
     os << '[' << NaNHelper(v.m_values[0]);
     for (size_t i = 1; i < N; ++i) {
@@ -105,7 +115,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr TinyVector operator+(const TinyVector& v) const
+  constexpr TinyVector
+  operator+(const TinyVector& v) const
   {
     TinyVector sum;
     for (size_t i = 0; i < N; ++i) {
@@ -115,7 +126,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr TinyVector operator+(TinyVector&& v) const
+  constexpr TinyVector
+  operator+(TinyVector&& v) const
   {
     for (size_t i = 0; i < N; ++i) {
       v.m_values[i] += m_values[i];
@@ -124,7 +136,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr TinyVector operator-(const TinyVector& v) const
+  constexpr TinyVector
+  operator-(const TinyVector& v) const
   {
     TinyVector difference;
     for (size_t i = 0; i < N; ++i) {
@@ -134,7 +147,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr TinyVector operator-(TinyVector&& v) const
+  constexpr TinyVector
+  operator-(TinyVector&& v) const
   {
     for (size_t i = 0; i < N; ++i) {
       v.m_values[i] = m_values[i] - v.m_values[i];
@@ -143,7 +157,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr TinyVector& operator+=(const TinyVector& v)
+  constexpr TinyVector&
+  operator+=(const TinyVector& v)
   {
     for (size_t i = 0; i < N; ++i) {
       m_values[i] += v.m_values[i];
@@ -152,15 +167,8 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr void operator+=(const volatile TinyVector& v) volatile
-  {
-    for (size_t i = 0; i < N; ++i) {
-      m_values[i] += v.m_values[i];
-    }
-  }
-
-  PUGS_INLINE
-  constexpr TinyVector& operator-=(const TinyVector& v)
+  constexpr TinyVector&
+  operator-=(const TinyVector& v)
   {
     for (size_t i = 0; i < N; ++i) {
       m_values[i] -= v.m_values[i];
@@ -169,21 +177,24 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr T& operator[](size_t i) noexcept(NO_ASSERT)
+  constexpr T&
+  operator[](size_t i) noexcept(NO_ASSERT)
   {
     Assert(i < N);
     return m_values[i];
   }
 
   PUGS_INLINE
-  constexpr const T& operator[](size_t i) const noexcept(NO_ASSERT)
+  constexpr const T&
+  operator[](size_t i) const noexcept(NO_ASSERT)
   {
     Assert(i < N);
     return m_values[i];
   }
 
   PUGS_INLINE
-  constexpr TinyVector& operator=(ZeroType) noexcept
+  constexpr TinyVector&
+  operator=(ZeroType) noexcept
   {
     static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
     for (size_t i = 0; i < N; ++i) {
@@ -231,7 +242,7 @@ class [[nodiscard]] TinyVector
   constexpr TinyVector(const TinyVector&) noexcept = default;
 
   PUGS_INLINE
-  constexpr TinyVector(TinyVector && v) noexcept = default;
+  constexpr TinyVector(TinyVector&& v) noexcept = default;
 
   PUGS_INLINE
   ~TinyVector() noexcept = default;
diff --git a/src/language/utils/ASTNodeDataType.cpp b/src/language/utils/ASTNodeDataType.cpp
index b89afd13493929a463f4585abb892e7937548906..6007b2b0641ec4ecb21005f9d86201b509a449be 100644
--- a/src/language/utils/ASTNodeDataType.cpp
+++ b/src/language/utils/ASTNodeDataType.cpp
@@ -163,32 +163,37 @@ getTupleDataType(const ASTNode& type_node)
 std::string
 dataTypeName(const ASTNodeDataType& data_type)
 {
-  std::string name;
   switch (data_type) {
-  case ASTNodeDataType::undefined_t:
-    name = "undefined";
-    break;
-  case ASTNodeDataType::bool_t:
-    name = "B";
-    break;
-  case ASTNodeDataType::unsigned_int_t:
-    name = "N";
-    break;
-  case ASTNodeDataType::int_t:
-    name = "Z";
-    break;
-  case ASTNodeDataType::double_t:
-    name = "R";
-    break;
-  case ASTNodeDataType::vector_t:
-    name = "R^" + stringify(data_type.dimension());
-    break;
-  case ASTNodeDataType::matrix_t:
-    name = "R^" + stringify(data_type.numberOfRows()) + "x" + stringify(data_type.numberOfColumns());
-    break;
-  case ASTNodeDataType::tuple_t:
-    name = "(" + dataTypeName(data_type.contentType()) + ")";
-    break;
+  case ASTNodeDataType::undefined_t: {
+    return "undefined";
+  }
+  case ASTNodeDataType::bool_t: {
+    return "B";
+  }
+  case ASTNodeDataType::unsigned_int_t: {
+    return "N";
+  }
+  case ASTNodeDataType::int_t: {
+    return "Z";
+  }
+  case ASTNodeDataType::double_t: {
+    return "R";
+  }
+  case ASTNodeDataType::vector_t: {
+    std::ostringstream data_type_name_vector;
+    data_type_name_vector << "R^" << data_type.dimension();
+    return data_type_name_vector.str();
+  }
+  case ASTNodeDataType::matrix_t: {
+    std::ostringstream data_type_name_matrix;
+    data_type_name_matrix << "R^" << data_type.numberOfRows() << "x" << data_type.numberOfColumns();
+    return data_type_name_matrix.str();
+  }
+  case ASTNodeDataType::tuple_t: {
+    std::ostringstream data_type_name_tuple;
+    data_type_name_tuple << "(" << dataTypeName(data_type.contentType()) << ")";
+    return data_type_name_tuple.str();
+  }
   case ASTNodeDataType::list_t: {
     std::ostringstream data_type_name_list;
     const auto& data_type_list = data_type.contentTypeList();
@@ -205,35 +210,38 @@ dataTypeName(const ASTNodeDataType& data_type)
           data_type_name_list << '*' << dataTypeName(*data_type_list[i]);
         }
       }
-      name = data_type_name_list.str();
+      return data_type_name_list.str();
     } else {
-      name = "void";
+      return "void";
     }
-    break;
-  }
-  case ASTNodeDataType::string_t:
-    name = "string";
-    break;
-  case ASTNodeDataType::typename_t:
-    name = dataTypeName(data_type.contentType());
-    break;
-  case ASTNodeDataType::type_name_id_t:
-    name = "type_name_id";
-    break;
-  case ASTNodeDataType::type_id_t:
-    name = data_type.nameOfTypeId();
-    break;
-  case ASTNodeDataType::function_t:
-    name = "function";
-    break;
-  case ASTNodeDataType::builtin_function_t:
-    name = "builtin_function";
-    break;
-  case ASTNodeDataType::void_t:
-    name = "void";
-    break;
-  }
-  return name;
+  }
+  case ASTNodeDataType::string_t: {
+    return "string";
+  }
+  case ASTNodeDataType::typename_t: {
+    return dataTypeName(data_type.contentType());
+  }
+  case ASTNodeDataType::type_name_id_t: {
+    return "type_name_id";
+  }
+  case ASTNodeDataType::type_id_t: {
+    return data_type.nameOfTypeId();
+  }
+  case ASTNodeDataType::function_t: {
+    return "function";
+  }
+  case ASTNodeDataType::builtin_function_t: {
+    return "builtin_function";
+  }
+  case ASTNodeDataType::void_t: {
+    return "void";
+  }
+  // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("unknown data type");
+  }
+    // LCOV_EXCL_STOP
+  }
 }
 
 std::string
diff --git a/src/mesh/ItemArrayUtils.hpp b/src/mesh/ItemArrayUtils.hpp
index abba5cfe2631bf6c472f3a4c60a5c0eec8703acb..b1c3ba8cdad569f73f91029c24fb1333e05037ee 100644
--- a/src/mesh/ItemArrayUtils.hpp
+++ b/src/mesh/ItemArrayUtils.hpp
@@ -53,7 +53,7 @@ min(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src < dst) {
         // cannot be reached if initial value is the min
@@ -152,7 +152,7 @@ max(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src > dst) {
         // cannot be reached if initial value is the max
@@ -248,7 +248,7 @@ sum(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       dst += src;
     }
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index 3bff4868c30677ef8d11db3a94cb493c5bcb07a0..9ef745b550e42add842ac1c822cf65b07e34dc2b 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -49,7 +49,7 @@ min(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src < dst) {
         // cannot be reached if initial value is the min
@@ -143,7 +143,7 @@ max(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src > dst) {
         // cannot be reached if initial value is the max
@@ -236,7 +236,7 @@ sum(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       dst += src;
     }
diff --git a/src/mesh/SubItemArrayPerItemUtils.hpp b/src/mesh/SubItemArrayPerItemUtils.hpp
index 1c5ad24683c5527e8028896348ba73ff8160af74..d2553590ab0f0dad219cdfa82491c05a70b5edbe 100644
--- a/src/mesh/SubItemArrayPerItemUtils.hpp
+++ b/src/mesh/SubItemArrayPerItemUtils.hpp
@@ -57,7 +57,7 @@ min(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_a
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src < dst) {
         // cannot be reached if initial array is the min
@@ -159,7 +159,7 @@ max(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_a
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src > dst) {
         // cannot be reached if initial array is the max
@@ -258,7 +258,7 @@ sum(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& item_array
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       dst += src;
     }
diff --git a/src/mesh/SubItemValuePerItemUtils.hpp b/src/mesh/SubItemValuePerItemUtils.hpp
index 8443047745680161747998a8cd114ffac8fe2647..304acae23a31e610695a2b309e72e7d8035d01a5 100644
--- a/src/mesh/SubItemValuePerItemUtils.hpp
+++ b/src/mesh/SubItemValuePerItemUtils.hpp
@@ -55,7 +55,7 @@ min(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_v
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src < dst) {
         // cannot be reached if initial value is the min
@@ -155,7 +155,7 @@ max(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_v
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src > dst) {
         // cannot be reached if initial value is the max
@@ -252,7 +252,7 @@ sum(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& item_value
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       dst += src;
     }
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 1fb2ce14cee3f42e6c27c262e3ba6deefb5b555d..6587a19f2c77d1f063f38765b4e768d92798f7a4 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -743,7 +743,7 @@ class DiscreteFunctionP0
 
     PUGS_INLINE
     void
-    join(volatile std::decay_t<data_type>& dst, const volatile std::decay_t<data_type>& src) const
+    join(std::decay_t<data_type>& dst, const std::decay_t<data_type>& src) const
     {
       dst += src;
     }
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 3ee64824a650e18714a9b4b8523afb78cc89f643..4943d4c08d983f67d1550262e34aa80030d49bf1 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -279,7 +279,7 @@ min(const Array<DataType>& array)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src < dst) {
         dst = src;
@@ -342,7 +342,7 @@ max(const Array<DataType>& array)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       if (src > dst) {
         dst = src;
@@ -402,7 +402,7 @@ sum(const Array<DataType>& array)
 
     PUGS_INLINE
     void
-    join(volatile data_type& dst, const volatile data_type& src) const
+    join(data_type& dst, const data_type& src) const
     {
       dst += src;
     }
diff --git a/tests/test_ASTNodeAffectationExpressionBuilder.cpp b/tests/test_ASTNodeAffectationExpressionBuilder.cpp
index f52a33c297bdea413b9f728498eecf8d67e96b30..39fcf659b208ac75ebbf198bddfe01b4c6ed06a1 100644
--- a/tests/test_ASTNodeAffectationExpressionBuilder.cpp
+++ b/tests/test_ASTNodeAffectationExpressionBuilder.cpp
@@ -173,6 +173,11 @@ const auto builtin_data_type = ast_node_data_type_from<std::shared_ptr<const dou
 
 #ifdef __clang__
 #pragma clang optimize off
+#else   // __clang__
+#ifdef __GNUG__
+#pragma GCC push_options
+#pragma GCC optimize("O0")
+#endif   // __GNUG__
 #endif   // __clang__
 
 // clazy:excludeall=non-pod-global-static
@@ -7530,4 +7535,8 @@ let a:void, a = 3;
 
 #ifdef __clang__
 #pragma clang optimize on
+#else   // __clang__
+#ifdef __GNUG__
+#pragma GCC pop_options
+#endif   // __GNUG__
 #endif   // __clang__
diff --git a/tests/test_ASTNodeFunctionExpressionBuilder.cpp b/tests/test_ASTNodeFunctionExpressionBuilder.cpp
index 7d760d2224fe4bbfb836b409a035830efd18b6fa..f26ee8ad39eac1c79150d3283ef71d45eba4ca71 100644
--- a/tests/test_ASTNodeFunctionExpressionBuilder.cpp
+++ b/tests/test_ASTNodeFunctionExpressionBuilder.cpp
@@ -111,6 +111,14 @@ const auto builtin_data_type = ast_node_data_type_from<std::shared_ptr<const dou
     REQUIRE_THROWS_WITH(ASTNodeExpressionBuilder{*ast}, error);                    \
   }
 
+#ifdef __clang__
+#else   // __clang__
+#ifdef __GNUG__
+#pragma GCC push_options
+#pragma GCC optimize("O0")
+#endif   // __GNUG__
+#endif   // __clang__
+
 // clazy:excludeall=non-pod-global-static
 
 TEST_CASE("ASTNodeFunctionExpressionBuilder", "[language]")
@@ -4065,3 +4073,10 @@ let h:R*void -> R, (x,void) -> 2.5;
     }
   }
 }
+
+#ifdef __clang__
+#else   // __clang__
+#ifdef __GNUG__
+#pragma GCC pop_options
+#endif   // __GNUG__
+#endif   // __clang__