diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index 48f94728a91c705d61fcfe537084ade59818f3f1..e64f2e1bdc0bc1e934aac22258461947c1c0ffd3 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -8,7 +8,7 @@
 #include <mesh/Synchronizer.hpp>
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/PugsTraits.hpp>
-#include <utils/ReproductibleSumMPI.hpp>
+#include <utils/ReproducibleSumMPI.hpp>
 
 #include <iostream>
 
@@ -294,20 +294,20 @@ sum(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
   auto is_owned = get_is_owned(*item_value.connectivity_ptr());
 
-  if (ReproductibleSumManager::reproductibleSums()) {
+  if (ReproducibleSumManager::reproducibleSums()) {
     if constexpr (std::is_floating_point_v<data_type>) {
-      ReproductibleScalarSum r_sum(item_value.arrayView(), is_owned.arrayView());
-      using ReproductibleSumType = decltype(r_sum);
-      return ReproductibleSumType::getValue(allReduceReproductibleSum(r_sum));
+      ReproducibleScalarSum r_sum(item_value.arrayView(), is_owned.arrayView());
+      using ReproducibleSumType = decltype(r_sum);
+      return ReproducibleSumType::getValue(allReduceReproducibleSum(r_sum));
     } else if constexpr (not std::is_arithmetic_v<data_type>) {
       if constexpr (is_tiny_vector_v<data_type> and std::is_floating_point_v<typename data_type::data_type>) {
-        ReproductibleTinyVectorSum r_sum(item_value.arrayView(), is_owned.arrayView());
-        using ReproductibleSumType = decltype(r_sum);
-        return ReproductibleSumType::getValue(allReduceReproductibleSum(r_sum));
+        ReproducibleTinyVectorSum r_sum(item_value.arrayView(), is_owned.arrayView());
+        using ReproducibleSumType = decltype(r_sum);
+        return ReproducibleSumType::getValue(allReduceReproducibleSum(r_sum));
       } else if constexpr (is_tiny_matrix_v<data_type> and std::is_floating_point_v<typename data_type::data_type>) {
-        ReproductibleTinyMatrixSum r_sum(item_value.arrayView(), is_owned.arrayView());
-        using ReproductibleSumType = decltype(r_sum);
-        return ReproductibleSumType::getValue(allReduceReproductibleSum(r_sum));
+        ReproducibleTinyMatrixSum r_sum(item_value.arrayView(), is_owned.arrayView());
+        using ReproducibleSumType = decltype(r_sum);
+        return ReproducibleSumType::getValue(allReduceReproducibleSum(r_sum));
       } else {
         const DataType local_sum = ItemValueSum{item_value, is_owned};
         return parallel::allReduceSum(local_sum);
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index fa3299bea5fefa99be94f4038ebae6da2c9779b8..2f0146e2b9d08c231eaadee7e7f50eeea439ddaf 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -6,8 +6,8 @@
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/PugsUtils.hpp>
-#include <utils/ReproductibleSumManager.hpp>
-#include <utils/ReproductibleSumUtils.hpp>
+#include <utils/ReproducibleSumManager.hpp>
+#include <utils/ReproducibleSumUtils.hpp>
 #include <utils/Types.hpp>
 
 #include <algorithm>
@@ -26,7 +26,8 @@ class [[nodiscard]] Array
     const size_t m_size;
 
    public:
-    friend std::ostream& operator<<(std::ostream& os, const UnsafeArrayView& x)
+    friend std::ostream&
+    operator<<(std::ostream& os, const UnsafeArrayView& x)
     {
       if (x.size() > 0) {
         os << 0 << ':' << NaNHelper(x[0]);
@@ -37,18 +38,21 @@ class [[nodiscard]] Array
       return os;
     }
 
-    [[nodiscard]] PUGS_INLINE size_t size() const
+    [[nodiscard]] PUGS_INLINE size_t
+    size() const
     {
       return m_size;
     }
 
-    [[nodiscard]] PUGS_INLINE DataType& operator[](size_t i) const
+    [[nodiscard]] PUGS_INLINE DataType&
+    operator[](size_t i) const
     {
       Assert(i < m_size, "invalid index");
       return m_values[i];
     }
 
-    PUGS_INLINE void fill(const DataType& data) const
+    PUGS_INLINE void
+    fill(const DataType& data) const
     {
       for (size_t i = 0; i < m_size; ++i) {
         m_values[i] = data;
@@ -56,7 +60,7 @@ class [[nodiscard]] Array
     }
 
     UnsafeArrayView& operator=(const UnsafeArrayView&) = delete;
-    UnsafeArrayView& operator=(UnsafeArrayView&&) = delete;
+    UnsafeArrayView& operator=(UnsafeArrayView&&)      = delete;
 
     UnsafeArrayView(const Array<DataType>& array, index_type begin, index_type size)
       : m_values{&array[begin]}, m_size{size}
@@ -80,12 +84,14 @@ class [[nodiscard]] Array
   friend Array<std::add_const_t<DataType>>;
 
  public:
-  [[nodiscard]] PUGS_INLINE size_t size() const noexcept
+  [[nodiscard]] PUGS_INLINE size_t
+  size() const noexcept
   {
     return m_values.extent(0);
   }
 
-  [[nodiscard]] friend PUGS_INLINE Array<std::remove_const_t<DataType>> copy(const Array<DataType>& source)
+  [[nodiscard]] friend PUGS_INLINE Array<std::remove_const_t<DataType>>
+  copy(const Array<DataType>& source)
   {
     Array<std::remove_const_t<DataType>> image(source.size());
     Kokkos::deep_copy(image.m_values, source.m_values);
@@ -93,8 +99,8 @@ class [[nodiscard]] Array
     return image;
   }
 
-  friend PUGS_INLINE void copy_to(const Array<DataType>& source,
-                                  const Array<std::remove_const_t<DataType>>& destination)
+  friend PUGS_INLINE void
+  copy_to(const Array<DataType>& source, const Array<std::remove_const_t<DataType>>& destination)
   {
     Assert(source.size() == destination.size(), "incompatible Array sizes");
     Kokkos::deep_copy(destination.m_values, source.m_values);
@@ -108,14 +114,16 @@ class [[nodiscard]] Array
                                                                       typename Array<DataType2>::index_type begin,
                                                                       typename Array<DataType2>::index_type size);
 
-  [[nodiscard]] PUGS_INLINE DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  [[nodiscard]] PUGS_INLINE DataType&
+  operator[](index_type i) const noexcept(NO_ASSERT)
   {
     Assert(i < m_values.extent(0), "invalid index");
     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 Array of const");
 
@@ -123,7 +131,8 @@ class [[nodiscard]] Array
   }
 
   template <typename DataType2>
-  PUGS_INLINE Array& operator=(const Array<DataType2>& array) noexcept
+  PUGS_INLINE Array&
+  operator=(const Array<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>>(),
@@ -166,7 +175,8 @@ class [[nodiscard]] Array
 #endif   // NDEBUG
   }
 
-  friend std::ostream& operator<<(std::ostream& os, const Array& x)
+  friend std::ostream&
+  operator<<(std::ostream& os, const Array& x)
   {
     if (x.size() > 0) {
       os << 0 << ':' << NaNHelper(x[0]);
@@ -184,13 +194,14 @@ class [[nodiscard]] Array
   Array(const Array&) = default;
 
   template <typename DataType2>
-  PUGS_INLINE Array(const Array<DataType2>& array) noexcept
+  PUGS_INLINE
+  Array(const Array<DataType2>& array) noexcept
   {
     this->operator=(array);
   }
 
   PUGS_INLINE
-  Array(Array &&) = default;
+  Array(Array&&) = default;
 
   PUGS_INLINE
   ~Array() = default;
@@ -421,16 +432,16 @@ sum(const Array<DataType>& array)
     ~ArraySum() = default;
   };
 
-  if (ReproductibleSumManager::reproductibleSums()) {
+  if (ReproducibleSumManager::reproducibleSums()) {
     if constexpr (std::is_floating_point_v<data_type>) {
-      ReproductibleScalarSum r_sum(array);
+      ReproducibleScalarSum r_sum(array);
       return r_sum.getSum();
     } else if constexpr (not std::is_arithmetic_v<data_type>) {
       if constexpr (is_tiny_vector_v<data_type> and std::is_floating_point_v<typename data_type::data_type>) {
-        ReproductibleTinyVectorSum r_sum(array);
+        ReproducibleTinyVectorSum r_sum(array);
         return r_sum.getSum();
       } else if constexpr (is_tiny_matrix_v<data_type> and std::is_floating_point_v<typename data_type::data_type>) {
-        ReproductibleTinyMatrixSum r_sum(array);
+        ReproducibleTinyMatrixSum r_sum(array);
         return r_sum.getSum();
       } else {
         return ArraySum(array);
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 3e61d14bb59cf96b873e13c55d92a6e6d9fe3f07..52daa45e8fc653c3cd45019b4dc213cbcd14cbb7 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -13,7 +13,7 @@ add_library(
   PETScWrapper.cpp
   PugsUtils.cpp
   RandomEngine.cpp
-  ReproductibleSumManager.cpp
+  ReproducibleSumManager.cpp
   RevisionInfo.cpp
   SignalManager.cpp
   SLEPcWrapper.cpp
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index 538d700745bb7f24749fa7e878c28a64db0e0bfd..9eeef02ad3f860be23154d963c2197204b3d9a8e 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -6,7 +6,7 @@
 #include <utils/FPEManager.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PETScWrapper.hpp>
-#include <utils/ReproductibleSumManager.hpp>
+#include <utils/ReproducibleSumManager.hpp>
 #include <utils/RevisionInfo.hpp>
 #include <utils/SLEPcWrapper.hpp>
 #include <utils/SignalManager.hpp>
@@ -118,8 +118,8 @@ initialize(int& argc, char* argv[])
     bool pause_on_error = false;
     app.add_flag("-p,--pause-on-error", pause_on_error, "Pause for debugging on unexpected error [default: false]");
 
-    bool reproductible_sums = true;
-    app.add_flag("--reproductible-sums,!--no-reproductible-sums", show_preamble,
+    bool reproducible_sums = true;
+    app.add_flag("--reproducible-sums,!--no-reproducible-sums", show_preamble,
                  "Special treatment of array sums to ensure reproducibility [default: true]");
 
     std::atexit([]() { std::cout << rang::style::reset; });
@@ -135,7 +135,7 @@ initialize(int& argc, char* argv[])
     ConsoleManager::setShowPreamble(show_preamble);
     ConsoleManager::init(enable_color);
     SignalManager::setPauseForDebug(pause_on_error);
-    ReproductibleSumManager::setReproductibleSums(reproductible_sums);
+    ReproducibleSumManager::setReproducibleSums(reproducible_sums);
   }
 
   PETScWrapper::initialize(argc, argv);
diff --git a/src/utils/ReproductibleSumMPI.hpp b/src/utils/ReproducibleSumMPI.hpp
similarity index 77%
rename from src/utils/ReproductibleSumMPI.hpp
rename to src/utils/ReproducibleSumMPI.hpp
index 91caafded09f11b31d29427d434c12704b5869e3..81bca2bd201b4b88f91b75439ccee9731ed5c707 100644
--- a/src/utils/ReproductibleSumMPI.hpp
+++ b/src/utils/ReproducibleSumMPI.hpp
@@ -1,14 +1,14 @@
-#ifndef REPRODUCTIBLE_SUM_MPI_HPP
-#define REPRODUCTIBLE_SUM_MPI_HPP
+#ifndef REPRODUCIBLE_SUM_MPI_HPP
+#define REPRODUCIBLE_SUM_MPI_HPP
 
 #include <utils/Messenger.hpp>
-#include <utils/ReproductibleSumUtils.hpp>
+#include <utils/ReproducibleSumUtils.hpp>
 #include <utils/pugs_config.hpp>
 
 #ifdef PUGS_HAS_MPI
 
 template <typename RSumType>
-struct ReproductibeMPIReducer
+struct ReproducibleMPIReducer
 {
   static void
   all_reduce_sum_embedder(void* local_value, void* sum_value, int* length, MPI_Datatype*)
@@ -24,7 +24,7 @@ struct ReproductibeMPIReducer
 
 template <typename ArrayT, typename BoolArrayT>
 auto
-allReduceReproductibleSum(const ReproductibleScalarSum<ArrayT, BoolArrayT>& s)
+allReduceReproducibleSum(const ReproducibleScalarSum<ArrayT, BoolArrayT>& s)
 {
   using DataType = std::decay_t<typename ArrayT::data_type>;
   using RSumType = std::decay_t<decltype(s)>;
@@ -38,7 +38,7 @@ allReduceReproductibleSum(const ReproductibleScalarSum<ArrayT, BoolArrayT>& s)
   MPI_Type_commit(&mpi_bin_type);
 
   MPI_Op mpi_bin_sum_op;
-  MPI_Op_create(ReproductibeMPIReducer<RSumType>::all_reduce_sum_embedder, 1, &mpi_bin_sum_op);
+  MPI_Op_create(ReproducibleMPIReducer<RSumType>::all_reduce_sum_embedder, 1, &mpi_bin_sum_op);
 
   BinType sum_bin = zero;
   MPI_Allreduce(&local_bin, &sum_bin, 1, mpi_bin_type, mpi_bin_sum_op, MPI_COMM_WORLD);
@@ -48,7 +48,7 @@ allReduceReproductibleSum(const ReproductibleScalarSum<ArrayT, BoolArrayT>& s)
 
 template <typename ArrayT, typename BoolArrayT>
 auto
-allReduceReproductibleSum(const ReproductibleTinyVectorSum<ArrayT, BoolArrayT>& s)
+allReduceReproducibleSum(const ReproducibleTinyVectorSum<ArrayT, BoolArrayT>& s)
 {
   using TinyVectorType = std::decay_t<typename ArrayT::data_type>;
   using DataType       = typename TinyVectorType::data_type;
@@ -63,7 +63,7 @@ allReduceReproductibleSum(const ReproductibleTinyVectorSum<ArrayT, BoolArrayT>&
   MPI_Type_commit(&mpi_bin_type);
 
   MPI_Op mpi_bin_sum_op;
-  MPI_Op_create(ReproductibeMPIReducer<RSumType>::all_reduce_sum_embedder, 1, &mpi_bin_sum_op);
+  MPI_Op_create(ReproducibleMPIReducer<RSumType>::all_reduce_sum_embedder, 1, &mpi_bin_sum_op);
 
   BinType sum_bin = zero;
   MPI_Allreduce(&local_bin, &sum_bin, 1, mpi_bin_type, mpi_bin_sum_op, MPI_COMM_WORLD);
@@ -73,7 +73,7 @@ allReduceReproductibleSum(const ReproductibleTinyVectorSum<ArrayT, BoolArrayT>&
 
 template <typename ArrayT, typename BoolArrayT>
 auto
-allReduceReproductibleSum(const ReproductibleTinyMatrixSum<ArrayT, BoolArrayT>& s)
+allReduceReproducibleSum(const ReproducibleTinyMatrixSum<ArrayT, BoolArrayT>& s)
 {
   using TinyVectorType = std::decay_t<typename ArrayT::data_type>;
   using DataType       = typename TinyVectorType::data_type;
@@ -88,7 +88,7 @@ allReduceReproductibleSum(const ReproductibleTinyMatrixSum<ArrayT, BoolArrayT>&
   MPI_Type_commit(&mpi_bin_type);
 
   MPI_Op mpi_bin_sum_op;
-  MPI_Op_create(ReproductibeMPIReducer<RSumType>::all_reduce_sum_embedder, 1, &mpi_bin_sum_op);
+  MPI_Op_create(ReproducibleMPIReducer<RSumType>::all_reduce_sum_embedder, 1, &mpi_bin_sum_op);
 
   BinType sum_bin = zero;
   MPI_Allreduce(&local_bin, &sum_bin, 1, mpi_bin_type, mpi_bin_sum_op, MPI_COMM_WORLD);
@@ -100,25 +100,25 @@ allReduceReproductibleSum(const ReproductibleTinyMatrixSum<ArrayT, BoolArrayT>&
 
 template <typename ArrayT, typename BoolArrayT>
 auto
-allReduceReproductibleSum(const ReproductibleScalarSum<ArrayT, BoolArrayT>& s)
+allReduceReproducibleSum(const ReproducibleScalarSum<ArrayT, BoolArrayT>& s)
 {
   return s.getSummationBin();
 }
 
 template <typename ArrayT, typename BoolArrayT>
 auto
-allReduceReproductibleSum(const ReproductibleTinyVectorSum<ArrayT, BoolArrayT>& s)
+allReduceReproducibleSum(const ReproducibleTinyVectorSum<ArrayT, BoolArrayT>& s)
 {
   return s.getSummationBin();
 }
 
 template <typename ArrayT, typename BoolArrayT>
 auto
-allReduceReproductibleSum(const ReproductibleTinyMatrixSum<ArrayT, BoolArrayT>& s)
+allReduceReproducibleSum(const ReproducibleTinyMatrixSum<ArrayT, BoolArrayT>& s)
 {
   return s.getSummationBin();
 }
 
 #endif   // PUGS_HAS_MPI
 
-#endif   // REPRODUCTIBLE_SUM_MPI_HPP
+#endif   // REPRODUCIBLE_SUM_MPI_HPP
diff --git a/src/utils/ReproducibleSumManager.cpp b/src/utils/ReproducibleSumManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b0f7f7d545e36346b7d0c3280024b81198cdce9c
--- /dev/null
+++ b/src/utils/ReproducibleSumManager.cpp
@@ -0,0 +1,15 @@
+#include <utils/ReproducibleSumManager.hpp>
+
+bool ReproducibleSumManager::s_reproducible_sums = true;
+
+bool
+ReproducibleSumManager::reproducibleSums()
+{
+  return s_reproducible_sums;
+}
+
+void
+ReproducibleSumManager::setReproducibleSums(bool reproducible_sum)
+{
+  s_reproducible_sums = reproducible_sum;
+}
diff --git a/src/utils/ReproducibleSumManager.hpp b/src/utils/ReproducibleSumManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..40e380fd6747a714906fbc15fc85667f551c8586
--- /dev/null
+++ b/src/utils/ReproducibleSumManager.hpp
@@ -0,0 +1,14 @@
+#ifndef REPRODUCIBLE_SUM_MANAGER_HPP
+#define REPRODUCIBLE_SUM_MANAGER_HPP
+
+class ReproducibleSumManager
+{
+ private:
+  static bool s_reproducible_sums;
+
+ public:
+  static void setReproducibleSums(bool);
+  static bool reproducibleSums();
+};
+
+#endif   // REPRODUCIBLE_SUM_MANAGER_HPP
diff --git a/src/utils/ReproductibleSumUtils.hpp b/src/utils/ReproducibleSumUtils.hpp
similarity index 88%
rename from src/utils/ReproductibleSumUtils.hpp
rename to src/utils/ReproducibleSumUtils.hpp
index e1f78a37fafc4522ffac80cec9047063dc3ddaf7..f2d52aa20122452e96130c0c29b6d12ccc32a4aa 100644
--- a/src/utils/ReproductibleSumUtils.hpp
+++ b/src/utils/ReproducibleSumUtils.hpp
@@ -1,5 +1,5 @@
-#ifndef REPRODUCTIBLE_SUM_UTILS_HPP
-#define REPRODUCTIBLE_SUM_UTILS_HPP
+#ifndef REPRODUCIBLE_SUM_UTILS_HPP
+#define REPRODUCIBLE_SUM_UTILS_HPP
 
 #include <utils/PugsUtils.hpp>
 #include <utils/Types.hpp>
@@ -7,7 +7,7 @@
 template <typename DataType>
 class Array;
 
-namespace reproductible_sum_utils
+namespace reproducible_sum_utils
 {
 template <size_t NumberOfBits>
 struct IntegerFromBitSize
@@ -102,24 +102,25 @@ constexpr inline size_t bin_max_size<float> = 1024;
 
 struct NoMask
 {
-  PUGS_INLINE bool operator[](size_t) const
+  PUGS_INLINE bool
+  operator[](size_t) const
   {
     return true;
   }
 };
 
-}   // namespace reproductible_sum_utils
+}   // namespace reproducible_sum_utils
 
-template <typename ArrayT, typename MaskType = reproductible_sum_utils::NoMask>
-class ReproductibleScalarSum
+template <typename ArrayT, typename MaskType = reproducible_sum_utils::NoMask>
+class ReproducibleScalarSum
 {
  public:
   using DataType = std::decay_t<typename ArrayT::data_type>;
   static_assert(std::is_floating_point_v<DataType>);
 
-  static constexpr size_t K     = reproductible_sum_utils::bin_number<DataType>;
-  static constexpr DataType eps = reproductible_sum_utils::bin_epsilon<DataType>;
-  static constexpr size_t W     = reproductible_sum_utils::bin_size<DataType>;
+  static constexpr size_t K     = reproducible_sum_utils::bin_number<DataType>;
+  static constexpr DataType eps = reproducible_sum_utils::bin_epsilon<DataType>;
+  static constexpr size_t W     = reproducible_sum_utils::bin_size<DataType>;
 
   struct Bin
   {
@@ -127,7 +128,7 @@ class ReproductibleScalarSum
     std::array<DataType, K> C;   // carry
 
     Bin& operator=(const Bin&) = default;
-    Bin& operator=(Bin&&) = default;
+    Bin& operator=(Bin&&)      = default;
 
     Bin(Bin&&)      = default;
     Bin(const Bin&) = default;
@@ -152,7 +153,7 @@ class ReproductibleScalarSum
   static void
   _shift(const size_t g, Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t k = K - 1; k >= g; --k) {
       bin.S[k] = bin.S[k - g];
@@ -169,7 +170,7 @@ class ReproductibleScalarSum
   {
     Assert(m > 0);
 
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
     if (m >= std::pow(DataType{2}, W - 1.) * ulp(bin.S[0])) {
       const size_t shift = 1 + std::floor(std::log2(m / (std::pow(DataType{2}, W - 1.) * ulp(bin.S[0]))) / W);
 
@@ -180,7 +181,7 @@ class ReproductibleScalarSum
   PUGS_INLINE
   void static _split2(DataType& S, DataType& x) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
     union
     {
       DataType as_DataType;
@@ -198,7 +199,7 @@ class ReproductibleScalarSum
   PUGS_INLINE static void
   _renormalize(Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t k = 0; k < K; ++k) {
       if (bin.S[k] >= 1.75 * ufp(bin.S[k])) {
@@ -218,7 +219,7 @@ class ReproductibleScalarSum
   static void
   addBinTo(Bin& bin, Bin& bin_sum)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     DataType ulp_bin = ulp(bin.S[0]);
     DataType ulp_sum = ulp(bin_sum.S[0]);
@@ -247,7 +248,7 @@ class ReproductibleScalarSum
   static DataType
   getValue(const Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     DataType value = 0;
     for (size_t k = 0; k < K; ++k) {
@@ -268,9 +269,9 @@ class ReproductibleScalarSum
     return getValue(m_summation_bin);
   }
 
-  ReproductibleScalarSum(const ArrayT& array, const MaskType mask = reproductible_sum_utils::NoMask{})
+  ReproducibleScalarSum(const ArrayT& array, const MaskType mask = reproducible_sum_utils::NoMask{})
   {
-    if constexpr (not std::is_same_v<MaskType, reproductible_sum_utils::NoMask>) {
+    if constexpr (not std::is_same_v<MaskType, reproducible_sum_utils::NoMask>) {
       static_assert(std::is_same_v<std::decay_t<typename MaskType::data_type>, bool>,
                     "when provided, mask must be an array of bool");
     }
@@ -278,7 +279,7 @@ class ReproductibleScalarSum
     using TeamPolicyT = Kokkos::TeamPolicy<Kokkos::IndexType<int>>;
     using TeamMemberT = TeamPolicyT::member_type;
 
-    int nx = reproductible_sum_utils::bin_max_size<DataType>;
+    int nx = reproducible_sum_utils::bin_max_size<DataType>;
     int ny = std::max(array.size() / nx, 1ul);
 
     const TeamPolicyT policy(ny, Kokkos::AUTO());
@@ -339,11 +340,11 @@ class ReproductibleScalarSum
     }
   }
 
-  ~ReproductibleScalarSum() = default;
+  ~ReproducibleScalarSum() = default;
 };
 
-template <typename ArrayT, typename MaskType = reproductible_sum_utils::NoMask>
-class ReproductibleTinyVectorSum
+template <typename ArrayT, typename MaskType = reproducible_sum_utils::NoMask>
+class ReproducibleTinyVectorSum
 {
  public:
   using TinyVectorType = std::decay_t<typename ArrayT::data_type>;
@@ -352,9 +353,9 @@ class ReproductibleTinyVectorSum
   using DataType = std::decay_t<typename TinyVectorType::data_type>;
   static_assert(std::is_floating_point_v<DataType>);
 
-  static constexpr size_t K     = reproductible_sum_utils::bin_number<DataType>;
-  static constexpr DataType eps = reproductible_sum_utils::bin_epsilon<DataType>;
-  static constexpr size_t W     = reproductible_sum_utils::bin_size<DataType>;
+  static constexpr size_t K     = reproducible_sum_utils::bin_number<DataType>;
+  static constexpr DataType eps = reproducible_sum_utils::bin_epsilon<DataType>;
+  static constexpr size_t W     = reproducible_sum_utils::bin_size<DataType>;
 
   struct Bin
   {
@@ -362,7 +363,7 @@ class ReproductibleTinyVectorSum
     std::array<TinyVectorType, K> C;   // carry
 
     Bin& operator=(const Bin&) = default;
-    Bin& operator=(Bin&&) = default;
+    Bin& operator=(Bin&&)      = default;
 
     Bin(Bin&&)      = default;
     Bin(const Bin&) = default;
@@ -390,7 +391,7 @@ class ReproductibleTinyVectorSum
   static void
   _shift(const size_t g, Bin& bin, const size_t& i_component) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t k = K - 1; k >= g; --k) {
       bin.S[k][i_component] = bin.S[k - g][i_component];
@@ -405,7 +406,7 @@ class ReproductibleTinyVectorSum
   PUGS_INLINE static void
   _update(const TinyVectorType& m, Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
     for (size_t i_component = 0; i_component < TinyVectorType::Dimension; ++i_component) {
       if (m[i_component] >= std::pow(DataType{2}, W - 1.) * ulp(bin.S[0][i_component])) {
         const size_t shift =
@@ -419,7 +420,7 @@ class ReproductibleTinyVectorSum
   PUGS_INLINE
   void static _split2(TinyVectorType& S, TinyVectorType& x) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
     union
     {
       DataType as_DataType;
@@ -442,7 +443,7 @@ class ReproductibleTinyVectorSum
   PUGS_INLINE static void
   _renormalize(Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t k = 0; k < K; ++k) {
       TinyVectorType& S_k = bin.S[k];
@@ -469,7 +470,7 @@ class ReproductibleTinyVectorSum
   static void
   addBinTo(Bin& bin, Bin& bin_sum)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t i_component = 0; i_component < TinyVectorType::Dimension; ++i_component) {
       DataType ulp_bin = ulp(bin.S[0][i_component]);
@@ -500,7 +501,7 @@ class ReproductibleTinyVectorSum
   static TinyVectorType
   getValue(const Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     TinyVectorType value = zero;
     for (size_t i_component = 0; i_component < TinyVectorType::Dimension; ++i_component) {
@@ -523,9 +524,9 @@ class ReproductibleTinyVectorSum
     return getValue(m_summation_bin);
   }
 
-  ReproductibleTinyVectorSum(const ArrayT& array, const MaskType mask = reproductible_sum_utils::NoMask{})
+  ReproducibleTinyVectorSum(const ArrayT& array, const MaskType mask = reproducible_sum_utils::NoMask{})
   {
-    if constexpr (not std::is_same_v<MaskType, reproductible_sum_utils::NoMask>) {
+    if constexpr (not std::is_same_v<MaskType, reproducible_sum_utils::NoMask>) {
       static_assert(std::is_same_v<std::decay_t<typename MaskType::data_type>, bool>,
                     "when provided, mask must be an array of bool");
     }
@@ -533,7 +534,7 @@ class ReproductibleTinyVectorSum
     using TeamPolicyT = Kokkos::TeamPolicy<Kokkos::IndexType<int>>;
     using TeamMemberT = TeamPolicyT::member_type;
 
-    int nx = reproductible_sum_utils::bin_max_size<DataType>;
+    int nx = reproducible_sum_utils::bin_max_size<DataType>;
     int ny = std::max(array.size() / nx, 1ul);
 
     const TeamPolicyT policy(ny, Kokkos::AUTO());
@@ -596,11 +597,11 @@ class ReproductibleTinyVectorSum
     }
   }
 
-  ~ReproductibleTinyVectorSum() = default;
+  ~ReproducibleTinyVectorSum() = default;
 };
 
-template <typename ArrayT, typename MaskType = reproductible_sum_utils::NoMask>
-class ReproductibleTinyMatrixSum
+template <typename ArrayT, typename MaskType = reproducible_sum_utils::NoMask>
+class ReproducibleTinyMatrixSum
 {
  public:
   using TinyMatrixType = std::decay_t<typename ArrayT::data_type>;
@@ -609,9 +610,9 @@ class ReproductibleTinyMatrixSum
   using DataType = std::decay_t<typename TinyMatrixType::data_type>;
   static_assert(std::is_floating_point_v<DataType>);
 
-  static constexpr size_t K     = reproductible_sum_utils::bin_number<DataType>;
-  static constexpr DataType eps = reproductible_sum_utils::bin_epsilon<DataType>;
-  static constexpr size_t W     = reproductible_sum_utils::bin_size<DataType>;
+  static constexpr size_t K     = reproducible_sum_utils::bin_number<DataType>;
+  static constexpr DataType eps = reproducible_sum_utils::bin_epsilon<DataType>;
+  static constexpr size_t W     = reproducible_sum_utils::bin_size<DataType>;
 
   struct Bin
   {
@@ -619,7 +620,7 @@ class ReproductibleTinyMatrixSum
     std::array<TinyMatrixType, K> C;   // carry
 
     Bin& operator=(const Bin&) = default;
-    Bin& operator=(Bin&&) = default;
+    Bin& operator=(Bin&&)      = default;
 
     Bin(Bin&&)      = default;
     Bin(const Bin&) = default;
@@ -648,7 +649,7 @@ class ReproductibleTinyMatrixSum
   static void
   _shift(const size_t g, Bin& bin, const size_t i_component, const size_t j_component) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t k = K - 1; k >= g; --k) {
       bin.S[k](i_component, j_component) = bin.S[k - g](i_component, j_component);
@@ -663,7 +664,7 @@ class ReproductibleTinyMatrixSum
   PUGS_INLINE static void
   _update(const TinyMatrixType& m, Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
     for (size_t i_component = 0; i_component < TinyMatrixType::NumberOfRows; ++i_component) {
       for (size_t j_component = 0; j_component < TinyMatrixType::NumberOfColumns; ++j_component) {
         if (m(i_component, j_component) >= std::pow(DataType{2}, W - 1.) * ulp(bin.S[0](i_component, j_component))) {
@@ -681,7 +682,7 @@ class ReproductibleTinyMatrixSum
   PUGS_INLINE
   void static _split2(TinyMatrixType& S, TinyMatrixType& x) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
     union
     {
       DataType as_DataType;
@@ -706,7 +707,7 @@ class ReproductibleTinyMatrixSum
   PUGS_INLINE static void
   _renormalize(Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t k = 0; k < K; ++k) {
       TinyMatrixType& S_k = bin.S[k];
@@ -735,7 +736,7 @@ class ReproductibleTinyMatrixSum
   static void
   addBinTo(Bin& bin, Bin& bin_sum)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     for (size_t i_component = 0; i_component < TinyMatrixType::NumberOfRows; ++i_component) {
       for (size_t j_component = 0; j_component < TinyMatrixType::NumberOfColumns; ++j_component) {
@@ -769,7 +770,7 @@ class ReproductibleTinyMatrixSum
   static TinyMatrixType
   getValue(const Bin& bin) noexcept(NO_ASSERT)
   {
-    using namespace reproductible_sum_utils;
+    using namespace reproducible_sum_utils;
 
     TinyMatrixType value = zero;
     for (size_t i_component = 0; i_component < TinyMatrixType::NumberOfRows; ++i_component) {
@@ -796,9 +797,9 @@ class ReproductibleTinyMatrixSum
     return getValue(m_summation_bin);
   }
 
-  ReproductibleTinyMatrixSum(const ArrayT& array, const MaskType mask = reproductible_sum_utils::NoMask{})
+  ReproducibleTinyMatrixSum(const ArrayT& array, const MaskType mask = reproducible_sum_utils::NoMask{})
   {
-    if constexpr (not std::is_same_v<MaskType, reproductible_sum_utils::NoMask>) {
+    if constexpr (not std::is_same_v<MaskType, reproducible_sum_utils::NoMask>) {
       static_assert(std::is_same_v<std::decay_t<typename MaskType::data_type>, bool>,
                     "when provided, mask must be an array of bool");
     }
@@ -806,7 +807,7 @@ class ReproductibleTinyMatrixSum
     using TeamPolicyT = Kokkos::TeamPolicy<Kokkos::IndexType<int>>;
     using TeamMemberT = TeamPolicyT::member_type;
 
-    int nx = reproductible_sum_utils::bin_max_size<DataType>;
+    int nx = reproducible_sum_utils::bin_max_size<DataType>;
     int ny = std::max(array.size() / nx, 1ul);
 
     const TeamPolicyT policy(ny, Kokkos::AUTO());
@@ -871,7 +872,7 @@ class ReproductibleTinyMatrixSum
     }
   }
 
-  ~ReproductibleTinyMatrixSum() = default;
+  ~ReproducibleTinyMatrixSum() = default;
 };
 
-#endif   // REPRODUCTIBLE_SUM_UTILS_HPP
+#endif   // REPRODUCIBLE_SUM_UTILS_HPP
diff --git a/src/utils/ReproductibleSumManager.cpp b/src/utils/ReproductibleSumManager.cpp
deleted file mode 100644
index 7f2851c5b740df50578d4f10800ad6255985d5d4..0000000000000000000000000000000000000000
--- a/src/utils/ReproductibleSumManager.cpp
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <utils/ReproductibleSumManager.hpp>
-
-bool ReproductibleSumManager::s_reproductible_sums = true;
-
-bool
-ReproductibleSumManager::reproductibleSums()
-{
-  return s_reproductible_sums;
-}
-
-void
-ReproductibleSumManager::setReproductibleSums(bool reproductible_sum)
-{
-  s_reproductible_sums = reproductible_sum;
-}
diff --git a/src/utils/ReproductibleSumManager.hpp b/src/utils/ReproductibleSumManager.hpp
deleted file mode 100644
index cd0b3a2df6006024ab4fa86e4ae88857c3a25233..0000000000000000000000000000000000000000
--- a/src/utils/ReproductibleSumManager.hpp
+++ /dev/null
@@ -1,14 +0,0 @@
-#ifndef REPRODUCTIBLE_SUM_MANAGER_HPP
-#define REPRODUCTIBLE_SUM_MANAGER_HPP
-
-class ReproductibleSumManager
-{
- private:
-  static bool s_reproductible_sums;
-
- public:
-  static void setReproductibleSums(bool);
-  static bool reproductibleSums();
-};
-
-#endif   // REPRODUCTIBLE_SUM_MANAGER_HPP
diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
index ff12fa0af4b211067d32fa71c9cfc85232aeb370..74dd2c37a96c84a7f72dc3d19138dae07d5ac6a2 100644
--- a/tests/test_Array.cpp
+++ b/tests/test_Array.cpp
@@ -14,7 +14,7 @@
 #include <valarray>
 #include <vector>
 
-#include <utils/ReproductibleSumUtils.hpp>
+#include <utils/ReproducibleSumUtils.hpp>
 #include <utils/SmallArray.hpp>
 #include <utils/Timer.hpp>
 
@@ -323,65 +323,131 @@ TEST_CASE("Array", "[utils]")
     }
   }
 
-  SECTION("reproducible double sum")
+  SECTION("reproducible floating point sum")
   {
-    Array<double> array(10'000);
+    SECTION("reproducible double sum")
+    {
+      Array<double> array(10'000);
 
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] = ((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = ((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+      }
+
+      const auto sum_before_shuffle = sum(array);
+
+      std::random_shuffle(&(array[0]), &(array[0]) + array.size());
+
+      const auto sum_after_shuffle = sum(array);
+
+      REQUIRE(sum_before_shuffle == sum_after_shuffle);
     }
 
-    const auto sum_before_shuffle = sum(array);
+    SECTION("reproducible float sum")
+    {
+      Array<float> array(10'000);
 
-    std::random_shuffle(&(array[0]), &(array[0]) + array.size());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = ((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+      }
 
-    const auto sum_after_shuffle = sum(array);
+      const auto sum_before_shuffle = sum(array);
 
-    REQUIRE(sum_before_shuffle == sum_after_shuffle);
-  }
+      std::random_shuffle(&(array[0]), &(array[0]) + array.size());
 
-  SECTION("reproducible TinyVector<3,float> sum")
-  {
-    Array<TinyVector<3, float>> array(10'000);
+      const auto sum_after_shuffle = sum(array);
 
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] =
-        TinyVector<3, float>(((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
-                             17 * i - 3,   //
-                             1);
+      REQUIRE(sum_before_shuffle == sum_after_shuffle);
     }
 
-    const auto sum_before_shuffle = sum(array);
+    SECTION("reproducible TinyVector<3,double> sum")
+    {
+      Array<TinyVector<3, double>> array(10'000);
 
-    std::random_shuffle(&(array[0]), &(array[0]) + array.size());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          TinyVector<3, double>(((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                17 * i - 3,   //
+                                1);
+      }
 
-    ReproductibleTinyVectorSum s0(array);
+      const auto sum_before_shuffle = sum(array);
 
-    const auto sum_after_shuffle = sum(array);
+      std::random_shuffle(&(array[0]), &(array[0]) + array.size());
 
-    REQUIRE(sum_before_shuffle == sum_after_shuffle);
-  }
+      ReproducibleTinyVectorSum s0(array);
 
-  SECTION("reproducible TinyMatrix<2, 3> sum")
-  {
-    Array<TinyMatrix<2, 3>> array(10'000);
+      const auto sum_after_shuffle = sum(array);
 
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] = TinyMatrix<2, 3>{((i < (array.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+      REQUIRE(sum_before_shuffle == sum_after_shuffle);
+    }
+
+    SECTION("reproducible TinyVector<3,float> sum")
+    {
+      Array<TinyVector<3, float>> array(10'000);
+
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          TinyVector<3, float>(((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                               17 * i - 3,   //
+                               1);
+      }
+
+      const auto sum_before_shuffle = sum(array);
+
+      std::random_shuffle(&(array[0]), &(array[0]) + array.size());
+
+      ReproducibleTinyVectorSum s0(array);
+
+      const auto sum_after_shuffle = sum(array);
+
+      REQUIRE(sum_before_shuffle == sum_after_shuffle);
+    }
+
+    SECTION("reproducible TinyMatrix<2, 3> sum")
+    {
+      Array<TinyMatrix<2, 3, double>> array(10'000);
+
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = TinyMatrix<2, 3, double>{((i < (array.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) *
+                                              std::sin(3 * i + 1),
+                                            17 * i - 3,            //
+                                            1,                     //
+                                            std::cos(7 * i) * i,   //
+                                            (2 - 3 * i) * 1E-30,   //
+                                            (i * 100 + 1000)};
+      }
+
+      const auto sum_before_shuffle = sum(array);
+
+      std::random_shuffle(&(array[0]), &(array[0]) + array.size());
+
+      const auto sum_after_shuffle = sum(array);
+
+      REQUIRE(sum_before_shuffle == sum_after_shuffle);
+    }
+
+    SECTION("reproducible TinyMatrix<2, 3, float> sum")
+    {
+      Array<TinyMatrix<2, 3, float>> array(10'000);
+
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          TinyMatrix<2, 3, float>(((i < (array.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
                                   17 * i - 3,            //
                                   1,                     //
                                   std::cos(7 * i) * i,   //
                                   (2 - 3 * i) * 1E-30,   //
-                                  (i * 100 + 1000)};
-    }
+                                  (i * 100 + 1000));
+      }
 
-    const auto sum_before_shuffle = sum(array);
+      const auto sum_before_shuffle = sum(array);
 
-    std::random_shuffle(&(array[0]), &(array[0]) + array.size());
+      std::random_shuffle(&(array[0]), &(array[0]) + array.size());
 
-    const auto sum_after_shuffle = sum(array);
+      const auto sum_after_shuffle = sum(array);
 
-    REQUIRE(sum_before_shuffle == sum_after_shuffle);
+      REQUIRE(sum_before_shuffle == sum_after_shuffle);
+    }
   }
 
   SECTION("checking for subArrayView")