diff --git a/src/main.cpp b/src/main.cpp
index 69c9ec388503ef93294a153f5a87e76bdccc8cb6..fbbdfd31e136c4318c9463382a245a6b714e1770 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,5 +1,5 @@
-#include <iostream>
 #include <PastisUtils.hpp>
+#include <PastisOStream.hpp>
 
 #include <rang.hpp>
 
@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
 
   try  {
     if (filename != "") {
-      std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
+      pout() << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
       Timer gmsh_timer;
       gmsh_timer.reset();
       GmshReader gmsh_reader(filename);
@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
                   break;
                 }
                 default: {
-                  std::cerr << "Unknown BCDescription\n";
+                  perr() << "Unknown BCDescription\n";
                   std::exit(1);
                 }
               }
@@ -134,7 +134,7 @@ int main(int argc, char *argv[])
           }
           vtk_writer.write(mesh, t, true); // forces last output
 
-          std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
+          pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
                     << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
           method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
@@ -195,7 +195,7 @@ int main(int argc, char *argv[])
                   break;
                 }
                 default: {
-                  std::cerr << "Unknown BCDescription\n";
+                  perr() << "Unknown BCDescription\n";
                   std::exit(1);
                 }
               }
@@ -241,7 +241,7 @@ int main(int argc, char *argv[])
           }
           vtk_writer.write(mesh, t, true); // forces last output
 
-          std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
+          pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
                     << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
           method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
@@ -291,7 +291,7 @@ int main(int argc, char *argv[])
                   break;
                 }
                 default: {
-                  std::cerr << "Unknown BCDescription\n";
+                  perr() << "Unknown BCDescription\n";
                   std::exit(1);
                 }
               }
@@ -336,7 +336,7 @@ int main(int argc, char *argv[])
           }
           vtk_writer.write(mesh, t, true); // forces last output
 
-          std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
+          pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
                     << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
           method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
@@ -344,15 +344,15 @@ int main(int argc, char *argv[])
         }
       }
 
-      std::cout << "* "  << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
+      pout() << "* "  << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
 
     } else {
-      std::cerr << "Connectivity1D defined by number of nodes no more implemented\n";
+      perr() << "Connectivity1D defined by number of nodes no more implemented\n";
       std::exit(0);
     }
   }
   catch (const AssertError& error) {
-    std::cerr << error << '\n';
+    perr() << error << '\n';
     std::exit(1);
   }
 
@@ -364,7 +364,7 @@ int main(int argc, char *argv[])
   }
 
   for (const auto& method_cost : method_cost_map) {
-    std::cout << "* ["
+    pout() << "* ["
               << rang::fgB::cyan
               << std::setw(size) << std::left
               << method_cost.first
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 90024775defe4ee51cbc6691328e4752ac1429b2..856a1452bf78ec323107142fc919b0bda4f79243 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -89,7 +89,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
         break;
       }
       default: {
-        std::cerr << "unexpected cell type!\n";
+        perr() << "unexpected cell type!\n";
         std::exit(0);
       }
     }
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index c4cd00ab7b8a3667379571f69a6281c8538be817..f876f70e11727c20fb72050da96a41bde8fe2f75 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -3,6 +3,8 @@
 
 #include <PastisMacros.hpp>
 #include <PastisAssert.hpp>
+
+#include <PastisOStream.hpp>
 #include <PastisUtils.hpp>
 
 #include <TinyVector.hpp>
@@ -136,7 +138,7 @@ class ConnectivityFace<3>
   friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
   {
     for (auto id : f.m_node_id_list) {
-      std::cout << id << ' ';
+      os << id << ' ';
     }
     return os;
   }
@@ -571,7 +573,7 @@ class Connectivity final
     const Face face(face_nodes);
     auto i_face = m_face_number_map.find(face);
     if (i_face == m_face_number_map.end()) {
-      std::cerr << "Face " << face << " not found!\n";
+      perr() << "Face " << face << " not found!\n";
       throw std::exception();
       std::exit(0);
     }
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 8c087e3baf65cec60ae10c5a130e51748dfe19e5..d6a77b13e6f124b825c54cbf422567281784fb79 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -16,13 +16,13 @@ computeConnectivityMatrix(const ConnectivityType& connectivity,
     const ConnectivityMatrix& child_to_item_matrix
         = connectivity._getMatrix(child_item_type, item_type);
 
-    std::cout << "computing connectivity "
+    pout() << "computing connectivity "
               << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
 
     item_to_child_item_matrix
         = this->_computeInverse(child_to_item_matrix);
   } else {
-    std::cerr << "unable to compute connectivity "
+    perr() << "unable to compute connectivity "
               << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
     std::exit(0);
   }
@@ -61,7 +61,7 @@ _computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
     size_t i=0;
     for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
       if (child_item_id != i) {
-        std::cerr << "sparse item numerotation NIY\n";
+        perr() << "sparse item numerotation NIY\n";
         std::exit(0);
       }
       ++i;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f187da556f28f8e4d0661077f6e397d8ce11b8f7..3a2d5c6ecc290a6ae71c234dd0551c7b48a22dca 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -96,30 +96,30 @@ void ErrorHandler::writeErrorMessage() const
 {
   switch(__type) {
     case asked: {
-      std::cerr << "\nremark: exit command explicitly called\n";
+      perr() << "\nremark: exit command explicitly called\n";
     }
     case normal: {
-      std::cerr << '\n' << __filename << ':' << __lineNumber
+      perr() << '\n' << __filename << ':' << __lineNumber
                 << ":remark: emitted the following message\n";
-      std::cerr << "error: " << __errorMessage << '\n';
+      perr() << "error: " << __errorMessage << '\n';
       break;
     }
     case compilation: {
-      std::cerr << "\nline " << __lineNumber << ':' << __errorMessage << '\n';
+      perr() << "\nline " << __lineNumber << ':' << __errorMessage << '\n';
       break;
     }
     case unexpected: {
-      std::cerr << '\n' << __filename << ':' << __lineNumber << ":\n" << __errorMessage << '\n';
-      std::cerr << "\nUNEXPECTED ERROR: this should not occure, please report it\n";
-      std::cerr << "\nBUG REPORT: Please send bug reports to:\n"
+      perr() << '\n' << __filename << ':' << __lineNumber << ":\n" << __errorMessage << '\n';
+      perr() << "\nUNEXPECTED ERROR: this should not occure, please report it\n";
+      perr() << "\nBUG REPORT: Please send bug reports to:\n"
                 << "  ff3d-dev@nongnu.org or freefem@ann.jussieu.fr\n"
                 << "or better, use the Bug Tracking System:\n"
                 << "  http://savannah.nongnu.org/bugs/?group=ff3d\n";
       break;
     }
     default: {
-      std::cerr << __filename << ':' << __lineNumber << ": " << __errorMessage << '\n';
-      std::cerr << __FILE__ << ':' << __LINE__ << ":remark:  error type not implemented!\n";
+      perr() << __filename << ':' << __lineNumber << ": " << __errorMessage << '\n';
+      perr() << __FILE__ << ':' << __LINE__ << ":remark:  error type not implemented!\n";
     }
   }
 }
@@ -144,7 +144,7 @@ GmshReader::GmshReader(const std::string& filename)
   try {
     m_fin.open(m_filename);
     if (not m_fin) {
-      std::cerr << "cannot read file '" << m_filename << "'\n";
+      perr() << "cannot read file '" << m_filename << "'\n";
       std::exit(0);
     }
 
@@ -206,7 +206,7 @@ GmshReader::GmshReader(const std::string& filename)
     __primitivesNames[14]    = "point";
     __supportedPrimitives[14]= true;
 
-    std::cout << "Reading file '" << m_filename << "'\n";
+    pout() << "Reading file '" << m_filename << "'\n";
 
     // Getting vertices list
     GmshReader::Keyword kw = this->__nextKeyword();
@@ -257,18 +257,18 @@ GmshReader::GmshReader(const std::string& filename)
           // 	}
           // }
 
-          //       std::cout << "- Binary format: ";
+          //       pout() << "- Binary format: ";
           // #ifdef    WORDS_BIGENDIAN
           //       if (not __convertEndian) {
-          // 	std::cout << "big endian\n";
+          // 	pout() << "big endian\n";
           //       } else {
-          // 	std::cout << "little endian\n";
+          // 	pout() << "little endian\n";
           //       }
           // #else  // WORDS_BIGENDIAN
           //       if (not __convertEndian) {
-          // 	std::cout << "little endian\n";
+          // 	pout() << "little endian\n";
           //       } else {
-          // 	std::cout << "big endian\n";
+          // 	pout() << "big endian\n";
           //       }
           // #endif // WORDS_BIGENDIAN
         }
@@ -304,7 +304,7 @@ GmshReader::GmshReader(const std::string& filename)
 void GmshReader::__readVertices()
 {
   const int numberOfVerices = this->_getInteger();
-  std::cout << "- Number Of Vertices: " << numberOfVerices << '\n';
+  pout() << "- Number Of Vertices: " << numberOfVerices << '\n';
   if (numberOfVerices<0) {
     throw ErrorHandler(__FILE__,__LINE__,
                        "reading file '"+this->m_filename
@@ -351,7 +351,7 @@ void GmshReader::__readVertices()
 // GmshReader::__readElements1()
 // {
 //   const int numberOfElements = this->_getInteger();
-//   std::cout << "- Number Of Elements: " << numberOfElements << '\n';
+//   pout() << "- Number Of Elements: " << numberOfElements << '\n';
 //   if (numberOfElements<0) {
 //     throw ErrorHandler(__FILE__,__LINE__,
 // 		       "reading file '"+m_filename
@@ -400,7 +400,7 @@ void
 GmshReader::__readElements2_2()
 {
   const int numberOfElements = this->_getInteger();
-  std::cout << "- Number Of Elements: " << numberOfElements << '\n';
+  pout() << "- Number Of Elements: " << numberOfElements << '\n';
   if (numberOfElements<0) {
     throw ErrorHandler(__FILE__,__LINE__,
                        "reading file '"+m_filename
@@ -516,7 +516,7 @@ __readPhysicalNames2_2()
     PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
     auto searched_physical_ref_id = m_physical_ref_map.find(physical_number);
     if (searched_physical_ref_id != m_physical_ref_map.end()) {
-      std::cerr << "Physical reference id '" << physical_ref_id
+      perr() << "Physical reference id '" << physical_ref_id
                 << "' already defined as '" << searched_physical_ref_id->second << "'!";
       std::exit(0);
     }
@@ -538,18 +538,18 @@ GmshReader::__proceedData()
 
   for (size_t i=0; i<elementNumber.dimension(); ++i) {
     if (elementNumber[i]>0) {
-      std::cout << "  - Number of "
+      pout() << "  - Number of "
                 << __primitivesNames[i]
                 << ": " << elementNumber[i];
       if (not(__supportedPrimitives[i])) {
-        std::cout << " [" << rang::fg::yellow << "IGNORED" << rang::style::reset << "]";
+        pout() << " [" << rang::fg::yellow << "IGNORED" << rang::style::reset << "]";
       } else {
-        std::cout << " references: ";
+        pout() << " references: ";
         for(std::set<size_t>::const_iterator iref
                 = elementReferences[i].begin();
             iref != elementReferences[i].end(); ++iref) {
-          if (iref != elementReferences[i].begin()) std::cout << ',';
-          std::cout << *iref;
+          if (iref != elementReferences[i].begin()) pout() << ',';
+          pout() << *iref;
         }
 
         switch (i) {
@@ -602,11 +602,11 @@ GmshReader::__proceedData()
           }
         }
       }
-      std::cout << '\n';
+      pout() << '\n';
     }
   }
 
-  std::cout << "- Building correspondance table\n";
+  pout() << "- Building correspondance table\n";
   int minNumber = std::numeric_limits<int>::max();
   int maxNumber = std::numeric_limits<int>::min();
   for (size_t i=0; i<__verticesNumbers.size(); ++i) {
@@ -785,10 +785,10 @@ GmshReader::__proceedData()
   dimension3_mask[12]=1;
   dimension3_mask[13]=1;
 
-  std::cout << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
-  std::cout << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
-  std::cout << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
-  std::cout << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
+  pout() << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
+  pout() << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
+  pout() << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
+  pout() << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
   if ((dimension3_mask, elementNumber)>0) {
     const size_t nb_cells = (dimension3_mask, elementNumber);
 
@@ -967,7 +967,7 @@ GmshReader::__proceedData()
     m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
 
   } else {
-    std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
+    perr() << "*** using a dimension 0 mesh is forbidden!\n";
     std::exit(0);
   }
 }
@@ -1005,7 +1005,7 @@ GmshReader::__nextKeyword()
 void GmshReader::
 __readGmshFormat2_2()
 {
-  std::cout << "- Reading Gmsh format 2.2\n";
+  pout() << "- Reading Gmsh format 2.2\n";
   GmshReader::Keyword kw = std::make_pair("", Unknown);
   while (kw.second != EndOfFile) {
     kw = this->__nextKeyword();
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 9f945e07b2c8b50889e014732c370cc74d1cb45a..0dab302789a506c9dfcaabd2e5ed3f4d8706cc2f 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -40,7 +40,7 @@ class MeshNodeBoundary
     parallel_for(face_list.size(), PASTIS_LAMBDA(const int& l){
         const auto& face_cells = face_to_cell_matrix[face_list[l]];
         if (face_cells.size()>1) {
-          std::cerr << "internal faces cannot be used to define mesh boundaries\n";
+          perr() << "internal faces cannot be used to define mesh boundaries\n";
           std::exit(1);
         }
       });
@@ -161,7 +161,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
   parallel_for(m_node_list.size(), PASTIS_LAMBDA(const size_t& r) {
       const R2& x = xr[m_node_list[r]];
       if ((x-origin,normal)>1E-13*length) {
-        std::cerr << "this FlatBoundary is not flat!\n";
+        perr() << "this FlatBoundary is not flat!\n";
         std::exit(1);
       }
     });
@@ -178,7 +178,7 @@ _getNormal(const MeshType& mesh)
   using R = TinyVector<1,double>;
 
   if (m_node_list.size() != 1) {
-    std::cerr << "Node boundaries in 1D require to have exactly 1 node\n";
+    perr() << "Node boundaries in 1D require to have exactly 1 node\n";
     std::exit(1);
   }
 
@@ -214,7 +214,7 @@ _getNormal(const MeshType& mesh)
   }
 
   if (xmin == xmax) {
-    std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
+    perr() << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
     std::exit(1);
   }
 
@@ -301,7 +301,7 @@ _getNormal(const MeshType& mesh)
 
 
   if (normal_l2 ==0) {
-    std::cerr << "not able to compute normal!\n";
+    perr() << "not able to compute normal!\n";
     std::exit(1);
   }
 
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 3afc434969159d4d8a547040be195669e2360f52..07ae8902dae5d00177e4f856e3fc5348db48bfc5 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -117,19 +117,19 @@ class AcousticSolver
     for (const auto& handler : m_boundary_condition_list) {
       switch (handler.boundaryCondition().type()) {
         case BoundaryCondition::normal_velocity: {
-          std::cerr << __FILE__ << ':' << __LINE__  << ": normal_velocity BC NIY\n";
+          perr() << __FILE__ << ':' << __LINE__  << ": normal_velocity BC NIY\n";
           std::exit(0);
           break;
         }
         case BoundaryCondition::velocity: {
-          std::cerr << __FILE__ << ':' << __LINE__  << ": velocity BC NIY\n";
+          perr() << __FILE__ << ':' << __LINE__  << ": velocity BC NIY\n";
           std::exit(0);
           break;
         }
         case BoundaryCondition::pressure: {
           // const PressureBoundaryCondition& pressure_bc
           //   = dynamic_cast<const PressureBoundaryCondition&>(handler.boundaryCondition());
-          std::cerr << __FILE__ << ':' << __LINE__ << ": pressure BC NIY\n";
+          perr() << __FILE__ << ':' << __LINE__ << ": pressure BC NIY\n";
           std::exit(0);
           break;
         }
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 1da880f33b95ff9a0696a7ef0ecab0f96860604a..5152bf5285ece836b5e8137755ab4935b434913e 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -9,6 +9,7 @@ add_library(
   BacktraceManager.cpp
   ConsoleManager.cpp
   FPEManager.cpp
+  PastisOStream.cpp
   PastisUtils.cpp
   RevisionInfo.cpp
   SignalManager.cpp)
diff --git a/src/utils/ConsoleManager.cpp b/src/utils/ConsoleManager.cpp
index b2a9984e49cef42cdb2b8b3464cc0d94d36e5744..cb1052d8e023b0b28f1df0986b2dca7cd6308fa9 100644
--- a/src/utils/ConsoleManager.cpp
+++ b/src/utils/ConsoleManager.cpp
@@ -1,4 +1,6 @@
 #include <ConsoleManager.hpp>
+#include <PastisOStream.hpp>
+
 #include <rang.hpp>
 
 bool ConsoleManager::isTerminal(std::ostream& os)
@@ -8,26 +10,26 @@ bool ConsoleManager::isTerminal(std::ostream& os)
 
 void ConsoleManager::init(const std::string& colorize)
 {
-  std::cout << "Console management: color ";
+  pout() << "Console management: color ";
   if (colorize == "auto") {
-    if (isTerminal(std::cout)) {
-      std::cout << rang::style::bold
+    if (isTerminal(pout())) {
+      pout() << rang::style::bold
 		<< rang::fgB::green
 		<< "enabled"
 		<< rang::fg::reset
 		<< rang::style::reset;
     } else {
-      std::cout << "disabled";
+      pout() << "disabled";
     }
-    std::cout << " [auto]\n";
+    pout() << " [auto]\n";
   } else if (colorize == "yes") {
     rang::setControlMode(rang::control::Force);
-    std::cout << rang::style::bold
+    pout() << rang::style::bold
 	      << rang::fgB::green
 	      << "enabled"
 	      << rang::fg::reset
 	      << rang::style::reset;
-    std::cout << " ["
+    pout() << " ["
 	      << rang::style::bold
 	      << rang::fgB::red
 	      << "forced"
@@ -36,9 +38,9 @@ void ConsoleManager::init(const std::string& colorize)
 	      << "]\n";
   } else if (colorize == "no") {
     rang::setControlMode(rang::control::Off);
-    std::cout << "disabled [forced]\n";
+    pout() << "disabled [forced]\n";
   } else {
-    std::cerr << "Unknown colorize option: " << colorize << '\n';
+    perr() << "Unknown colorize option: " << colorize << '\n';
     std::exit(1);
   }
 }
diff --git a/src/utils/FPEManager.cpp b/src/utils/FPEManager.cpp
index d4a4247682afa9f72c771933dcf091e00010d4bc..d852267132eda5917feaeba981ecbb2fdd5af1b9 100644
--- a/src/utils/FPEManager.cpp
+++ b/src/utils/FPEManager.cpp
@@ -1,5 +1,7 @@
 #include <FPEManager.hpp>
 #include <PastisMacros.hpp>
+#include <PastisOStream.hpp>
+
 #include <pastis_config.hpp>
 #include <rang.hpp>
 
@@ -57,7 +59,7 @@ int fedisableexcept(unsigned int excepts)
 
 void FPEManager::enable()
 {
-  std::cout << "FE management: "
+  pout() << "FE management: "
 	    << rang::style::bold
 	    << rang::fgB::green
 	    << "enabled"
@@ -68,7 +70,7 @@ void FPEManager::enable()
 
 void FPEManager::disable()
 {
-  std::cout << "FE management: "
+  pout() << "FE management: "
 	    << rang::style::bold
 	    << rang::fgB::red
 	    << "disabled"
@@ -81,7 +83,7 @@ void FPEManager::disable()
 
 void FPEManager::enable()
 {
-  std::cout << "FE management: enabled "
+  pout() << "FE management: enabled "
 	    << rang::fg::red
 	    << "[not supported]"
 	    << rang::fg::reset << '\n';
@@ -89,7 +91,7 @@ void FPEManager::enable()
 
 void FPEManager::disable()
 {
-  std::cout << "FE management: disable "
+  pout() << "FE management: disable "
 	    << rang::fg::red
 	    << "[not supported]"
 	    << rang::fg::reset << '\n';
diff --git a/src/utils/PastisOStream.cpp b/src/utils/PastisOStream.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7276b868ccdfa7faef9b27f231585a0883167b6
--- /dev/null
+++ b/src/utils/PastisOStream.cpp
@@ -0,0 +1,4 @@
+#include <PastisOStream.hpp>
+
+PastisOStream pout(std::cout);
+PastisOStream perr(std::cerr);
diff --git a/src/utils/PastisOStream.hpp b/src/utils/PastisOStream.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d8ec83072a870776fa23e363058fa42f0241151
--- /dev/null
+++ b/src/utils/PastisOStream.hpp
@@ -0,0 +1,41 @@
+#ifndef PASTIS_OSTREAM_HPP
+#define PASTIS_OSTREAM_HPP
+
+#include <PastisMacros.hpp>
+#include <iostream>
+
+class PastisOStream
+{
+ private:
+  std::ostream* m_os = nullptr;
+
+ public:
+  PASTIS_INLINE
+  std::ostream& operator()()
+  {
+    return *m_os;
+  }
+
+  void setOutput(std::ostream& os)
+  {
+    m_os = &os;
+  }
+
+  PastisOStream(std::ostream& os)
+      : m_os(&os)
+  {
+    ;
+  }
+
+  PastisOStream& operator=(const PastisOStream&) = delete;
+  PastisOStream& operator=(PastisOStream&&) = delete;
+  PastisOStream(const PastisOStream&) = delete;
+  PastisOStream(PastisOStream&&) = delete;
+
+  ~PastisOStream() = default;
+};
+
+extern PastisOStream pout;
+extern PastisOStream perr;
+
+#endif // PASTIS_OSTREAM_HPP
diff --git a/src/utils/PastisUtils.cpp b/src/utils/PastisUtils.cpp
index 4fea1e5d8e2ff8adf6d144123558601d4e67b29c..b7250c194bd17c1a9915cd41323a7b3dd587fc11 100644
--- a/src/utils/PastisUtils.cpp
+++ b/src/utils/PastisUtils.cpp
@@ -1,4 +1,6 @@
 #include <PastisUtils.hpp>
+#include <PastisOStream.hpp>
+
 #include <Kokkos_Core.hpp>
 
 #include <RevisionInfo.hpp>
@@ -17,35 +19,35 @@ std::string initialize(int& argc, char* argv[])
   long unsigned number = 10;
   std::string filename;
 
-  std::cout << "Pastis version: "
+  pout() << "Pastis version: "
             << rang::style::bold << RevisionInfo::version() << rang::style::reset << '\n';
 
-  std::cout << "-------------------- "
+  pout() << "-------------------- "
             << rang::fg::green
             << "git info"
             << rang::fg::reset
             <<" -------------------------"
             << '\n';
-  std::cout << "tag:  " << rang::style::bold << RevisionInfo::gitTag() << rang::style::reset << '\n';
-  std::cout << "HEAD: " << rang::style::bold << RevisionInfo::gitHead() << rang::style::reset << '\n';
-  std::cout << "hash: " << rang::style::bold << RevisionInfo::gitHash() << rang::style::reset << "  (";
+  pout() << "tag:  " << rang::style::bold << RevisionInfo::gitTag() << rang::style::reset << '\n';
+  pout() << "HEAD: " << rang::style::bold << RevisionInfo::gitHead() << rang::style::reset << '\n';
+  pout() << "hash: " << rang::style::bold << RevisionInfo::gitHash() << rang::style::reset << "  (";
 
   if (RevisionInfo::gitIsClean()) {
-    std::cout << rang::fgB::green << "clean" << rang::fg::reset;
+    pout() << rang::fgB::green << "clean" << rang::fg::reset;
   } else {
-    std::cout << rang::fgB::red << "dirty" << rang::fg::reset;
+    pout() << rang::fgB::red << "dirty" << rang::fg::reset;
   }
-  std::cout << ")\n";
-  std::cout << "-------------------- "
+  pout() << ")\n";
+  pout() << "-------------------- "
             << rang::fg::green
             << "build info"
             << rang::fg::reset
             <<" -----------------------"
             << '\n';
-  std::cout << "type:     " << rang::style::bold << BuildInfo::type() << rang::style::reset << '\n';
-  std::cout << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n';
-  std::cout << "devices:  " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
-  std::cout << "-------------------------------------------------------\n";
+  pout() << "type:     " << rang::style::bold << BuildInfo::type() << rang::style::reset << '\n';
+  pout() << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n';
+  pout() << "devices:  " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
+  pout() << "-------------------------------------------------------\n";
   {
     CLI::App app{"Pastis help"};
 
@@ -67,7 +69,7 @@ std::string initialize(int& argc, char* argv[])
     std::string pause_on_error="auto";
     app.add_set("--pause-on-error", pause_on_error, {"auto", "yes", "no"}, "Pause for debugging on unexpected error", true);
 
-    std::atexit([](){std::cout << rang::style::reset;});
+    std::atexit([](){pout() << rang::style::reset;});
     try {
       app.parse(argc, argv);
     } catch (const CLI::ParseError &e) {
@@ -81,17 +83,17 @@ std::string initialize(int& argc, char* argv[])
   }
 
   Kokkos::initialize(argc,argv);
-  std::cout << "-------------------- "
+  pout() << "-------------------- "
             << rang::fg::green
             << "exec info"
             << rang::fg::reset
             <<" ------------------------"
             << '\n';
 
-  std::cout << rang::style::bold;
-  Kokkos::DefaultExecutionSpace::print_configuration(std::cout);
-  std::cout << rang::style::reset;
-  std::cout << "-------------------------------------------------------\n";
+  pout() << rang::style::bold;
+  Kokkos::DefaultExecutionSpace::print_configuration(pout());
+  pout() << rang::style::reset;
+  pout() << "-------------------------------------------------------\n";
 
   return filename;
 }
diff --git a/src/utils/SignalManager.cpp b/src/utils/SignalManager.cpp
index 3e41a12ad47e0b5229cef9f1c7acc4fd3fa19822..09a1c5ab8f9a6b56f58e61cd645ba3da86db8a5d 100644
--- a/src/utils/SignalManager.cpp
+++ b/src/utils/SignalManager.cpp
@@ -1,5 +1,7 @@
 #include <SignalManager.hpp>
 
+#include <PastisOStream.hpp>
+
 #include <BacktraceManager.hpp>
 #include <ConsoleManager.hpp>
 
@@ -31,20 +33,20 @@ std::string SignalManager::signalName(const int& signal)
 void SignalManager::pauseForDebug(const int& signal)
 {
 
-  std::cout  << "\n======================================\n";
-  std::cout << rang::style::reset
-	    << rang::fg::reset
-	    << rang::style::bold;
-  std::cout << "to attach gdb to this process run\n";
-  std::cout << "\tgdb -pid "
-	    << rang::fg::red
-	    << getpid()
-	    << rang::fg::reset
-	    << '\n';
-  std::cout << "else press Control-C to exit\n";
-  std::cout << rang::style::reset;
-  std::cout  << "======================================\n";
-  std::cout << std::flush;
+  pout() << "\n======================================\n";
+  pout() << rang::style::reset
+       << rang::fg::reset
+       << rang::style::bold;
+  pout() << "to attach gdb to this process run\n";
+  pout() << "\tgdb -pid "
+       << rang::fg::red
+       << getpid()
+       << rang::fg::reset
+       << '\n';
+  pout() << "else press Control-C to exit\n";
+  pout() << rang::style::reset;
+  pout() << "======================================\n";
+  pout() << std::flush;
 
   pause();
   std::exit(signal);
@@ -58,11 +60,11 @@ void SignalManager::handler(int signal)
   std::signal(SIGINT,  SIG_DFL);
   std::signal(SIGABRT, SIG_DFL);
 
-  std::cerr << "\n *** "
+  perr() << "\n *** "
 	    << rang::style::reset
 	    << rang::fg::reset
 	    << rang::style::bold;
-  std::cerr << "Signal "
+  perr() << "Signal "
     	    << rang::fgB::red
 	    << signalName(signal)
     	    << rang::fg::reset
@@ -71,9 +73,9 @@ void SignalManager::handler(int signal)
 	    << " ***\n";
 
   BacktraceManager bm;
-  std::cerr << bm << '\n';
+  perr() << bm << '\n';
 
-  if ((ConsoleManager::isTerminal(std::cout) and
+  if ((ConsoleManager::isTerminal(pout()) and
        (s_pause_on_error == "auto")) or
       (s_pause_on_error == "yes")) {
     SignalManager::pauseForDebug(signal);
@@ -93,14 +95,14 @@ void SignalManager::init(const bool& enable)
     std::signal(SIGABRT, SignalManager::handler);
     std::signal(SIGPIPE, SignalManager::handler);
 
-    std::cout << "Signal management: "
+    pout() << "Signal management: "
 	      << rang::style::bold
 	      << rang::fgB::green
 	      << "enabled"
 	      << rang::fg::reset
 	      << rang::style::reset << '\n';
   } else {
-    std::cout << "Signal management: "
+    pout() << "Signal management: "
 	      << rang::style::bold
 	      << rang::fgB::red
 	      << "disabled"