diff --git a/.gitignore b/.gitignore
index cd78c3d5b..28bc43723 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,11 @@ __pycache__
 .cache
 *.old
 
+# Editor swap files
+*.swp
+*.swo
+*.swn
+
 # Created by https://www.toptal.com/developers/gitignore/api/cmake,linux,macos,visualstudiocode,python
 # Edit at https://www.toptal.com/developers/gitignore?templates=cmake,linux,macos,visualstudiocode,python
 
diff --git a/demos/FiniteVolume/AMR_Burgers_Hat.cpp b/demos/FiniteVolume/AMR_Burgers_Hat.cpp
index 0a975c4e8..fa514aa58 100644
--- a/demos/FiniteVolume/AMR_Burgers_Hat.cpp
+++ b/demos/FiniteVolume/AMR_Burgers_Hat.cpp
@@ -12,6 +12,7 @@
 #include <samurai/field.hpp>
 #include <samurai/io/hdf5.hpp>
 #include <samurai/io/restart.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "stencil_field.hpp"
@@ -269,7 +270,7 @@ int main(int argc, char* argv[])
         std::size_t ite_adapt = 0;
         while (true)
         {
-            std::cout << "\tmesh adaptation: " << ite_adapt++ << std::endl;
+            samurai::io::print("mesh adaptation: {}\n", ite_adapt++);
             samurai::update_ghost(phi);
             tag.resize();
             AMR_criteria(phi, tag);
@@ -287,7 +288,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         // Numerical scheme
         samurai::update_ghost(phi);
diff --git a/demos/FiniteVolume/BZ/bz_2d.cpp b/demos/FiniteVolume/BZ/bz_2d.cpp
index 05d7efd64..55c193729 100644
--- a/demos/FiniteVolume/BZ/bz_2d.cpp
+++ b/demos/FiniteVolume/BZ/bz_2d.cpp
@@ -9,6 +9,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/stencil_field.hpp>
 
 #include "../../LBM/boundary_conditions.hpp"
@@ -313,7 +314,7 @@ int main()
 
     while (t < Tf)
     {
-        fmt::print(fmt::format("Iteration = {:4d}, t: {}\n", nb_ite, t));
+        samurai::io::print(samurai::io::root, "{}", fmt::format("Iteration = {:4d}, t: {}\n", nb_ite, t));
 
         if (max_level > min_level)
         {
@@ -323,12 +324,12 @@ int main()
         tic();
         reaction(field, t, t + .5 * dt);
         auto duration = toc();
-        fmt::print(fmt::format("first reaction: {}\n", duration));
+        samurai::io::print("{}", fmt::format("first reaction: {}\n", duration));
 
         tic();
         RK4(field, dt, std::ceil(dt / dt_diffusion), update_bc, D_b, D_c);
         duration = toc();
-        fmt::print(fmt::format("diffusion: {}\n", duration));
+        samurai::io::print("{}", fmt::format("diffusion: {}\n", duration));
         /*
         rock4_integration(double tini, double tend, int neq, double *u,
                     func_rock fcn, double tol, int *info)
@@ -336,7 +337,7 @@ int main()
         tic();
         reaction(field, t + .5 * dt, t + dt);
         duration = toc();
-        fmt::print(fmt::format("second reaction: {}\n", duration));
+        samurai::io::print("{}", fmt::format("second reaction: {}\n", duration));
 
         if (nsave == 20)
         {
diff --git a/demos/FiniteVolume/BZ/bz_2d_AMR.cpp b/demos/FiniteVolume/BZ/bz_2d_AMR.cpp
index 4cbde0f21..700e4e907 100644
--- a/demos/FiniteVolume/BZ/bz_2d_AMR.cpp
+++ b/demos/FiniteVolume/BZ/bz_2d_AMR.cpp
@@ -10,6 +10,7 @@
 #include <samurai/mr/operators.hpp>
 #include <samurai/numeric/prediction.hpp>
 #include <samurai/numeric/projection.hpp>
+#include <samurai/print.hpp>
 #include <samurai/static_algorithm.hpp>
 #include <samurai/stencil_field.hpp>
 
@@ -691,7 +692,7 @@ int main()
 
     while (t < Tf)
     {
-        fmt::print(fmt::format("Iteration = {:4d}, t: {}\n", nb_ite, t));
+        samurai::io::print(samurai::io::root, "{}", fmt::format("Iteration = {:4d}, t: {}\n", nb_ite, t));
 
         std::size_t idx = 0;
 
@@ -710,12 +711,12 @@ int main()
         tic();
         reaction(field, t, t + .5 * dt);
         auto duration = toc();
-        fmt::print(fmt::format("first reaction: {}\n", duration));
+        samurai::io::print("{}", fmt::format("first reaction: {}\n", duration));
 
         tic();
         RK4(field, dt, std::ceil(dt / dt_diffusion), update_bc_for_level, D_b, D_c);
         duration = toc();
-        fmt::print(fmt::format("diffusion: {}\n", duration));
+        samurai::io::print("{}", fmt::format("diffusion: {}\n", duration));
         /*
         rock4_integration(double tini, double tend, int neq, double *u,
                     func_rock fcn, double tol, int *info)
@@ -723,7 +724,7 @@ int main()
         tic();
         reaction(field, t + .5 * dt, t + dt);
         duration = toc();
-        fmt::print(fmt::format("second reaction: {}\n", duration));
+        samurai::io::print("{}", fmt::format("second reaction: {}\n", duration));
 
         if (nsave == 20)
         {
diff --git a/demos/FiniteVolume/advection_1d.cpp b/demos/FiniteVolume/advection_1d.cpp
index b2601d2fb..3f3ae82ed 100644
--- a/demos/FiniteVolume/advection_1d.cpp
+++ b/demos/FiniteVolume/advection_1d.cpp
@@ -145,7 +145,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         samurai::update_ghost_mr(u);
         unp1.resize();
diff --git a/demos/FiniteVolume/advection_2d.cpp b/demos/FiniteVolume/advection_2d.cpp
index 89beea56f..c7f51db96 100644
--- a/demos/FiniteVolume/advection_2d.cpp
+++ b/demos/FiniteVolume/advection_2d.cpp
@@ -139,7 +139,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         samurai::update_ghost_mr(u);
         unp1.resize();
diff --git a/demos/FiniteVolume/advection_2d_user_bc.cpp b/demos/FiniteVolume/advection_2d_user_bc.cpp
index b8ed6f850..eb860634b 100644
--- a/demos/FiniteVolume/advection_2d_user_bc.cpp
+++ b/demos/FiniteVolume/advection_2d_user_bc.cpp
@@ -171,7 +171,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         samurai::update_ghost_mr(u);
         unp1.resize();
diff --git a/demos/FiniteVolume/burgers.cpp b/demos/FiniteVolume/burgers.cpp
index c965fe9eb..8854ab75d 100644
--- a/demos/FiniteVolume/burgers.cpp
+++ b/demos/FiniteVolume/burgers.cpp
@@ -56,7 +56,7 @@ int main_dim(int argc, char* argv[])
     using Box     = samurai::Box<double, dim>;
     using point_t = typename Box::point_t;
 
-    std::cout << "------------------------- Burgers -------------------------" << std::endl;
+    samurai::io::print("------------------------- Burgers -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -177,7 +177,7 @@ int main_dim(int argc, char* argv[])
         }
         else
         {
-            std::cerr << "Unmanaged initial solution '" << init_sol << "'.";
+            samurai::io::eprint("Unmanaged initial solution '{}' .\n", init_sol);
             return EXIT_FAILURE;
         }
     }
@@ -245,7 +245,7 @@ int main_dim(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -287,7 +287,7 @@ int main_dim(int argc, char* argv[])
                                              {
                                                  return exact_solution<n_comp>(coord, t);
                                              });
-            std::cout << ", L2-error: " << std::scientific << std::setprecision(2) << error;
+            samurai::io::print(", L2-error: {:.2e}", error);
 
             if (mesh.min_level() != mesh.max_level())
             {
@@ -298,19 +298,18 @@ int main_dim(int argc, char* argv[])
                                           {
                                               return exact_solution<n_comp>(coord, t);
                                           });
-                std::cout << ", L2-error (recons): " << std::scientific << std::setprecision(2) << error;
+                samurai::io::print(", L2-error (recons): {:.2e}", error);
             }
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     if constexpr (dim == 1)
     {
-        std::cout << std::endl;
-        std::cout << "Run the following command to view the results:" << std::endl;
-        std::cout << "python <<path to samurai>>/python/read_mesh.py " << filename << "_ite_ --field u level --start 0 --end " << nsave
-                  << std::endl;
+        samurai::io::print("\n");
+        samurai::io::print("Run the following command to view the results:\n");
+        samurai::io::print("python <<path to samurai>>/python/read_mesh.py {}_ite_ --field u level --start 0 --end {}\n", filename, nsave);
     }
 
     samurai::finalize();
diff --git a/demos/FiniteVolume/burgers_mra.cpp b/demos/FiniteVolume/burgers_mra.cpp
index ea6e74dbe..25957315e 100644
--- a/demos/FiniteVolume/burgers_mra.cpp
+++ b/demos/FiniteVolume/burgers_mra.cpp
@@ -82,7 +82,7 @@ void run_simulation(Field& u,
 {
     auto& mesh = u.mesh();
 
-    std::cout << std::endl << "max-level-flux enabled: " << scheme.enable_max_level_flux() << std::endl;
+    samurai::io::print("\nmax-level-flux enabled: {}\n", scheme.enable_max_level_flux());
 
     if (scheme.enable_max_level_flux())
     {
@@ -162,7 +162,7 @@ void run_simulation(Field& u,
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {:.2f}, dt = {}\n", nt++, t, dt);
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -179,8 +179,7 @@ void run_simulation(Field& u,
         auto u_recons = samurai::reconstruction(u);
 
         // Error
-        std::cout << ", L2-error: " << std::scientific;
-        std::cout.precision(2);
+        samurai::io::print(", L2-error: ");
         if (init_sol == "hat")
         {
             double error = samurai::L2_error(u,
@@ -188,14 +187,14 @@ void run_simulation(Field& u,
                                              {
                                                  return hat_exact_solution(coord[0], t);
                                              });
-            std::cout << "[w.r.t. exact, no recons] " << error;
+            samurai::io::print("[w.r.t. exact, no recons] {}", error);
 
             error = samurai::L2_error(u_recons,
                                       [&](const auto& coord)
                                       {
                                           return hat_exact_solution(coord[0], t);
                                       });
-            std::cout << ", [w.r.t. exact, recons] " << error;
+            samurai::io::print(", [w.r.t. exact, recons] {}", error);
         }
 
         double error = 0;
@@ -206,7 +205,7 @@ void run_simulation(Field& u,
                                    error += std::pow(u_max[cell] - u_recons[cell2], 2) * std::pow(cell.length, 1);
                                });
         error = std::sqrt(error);
-        std::cout << ", [w.r.t. max level, recons] " << error;
+        samurai::io::print(", [w.r.t. max level, recons] {}", error);
 
         // Save the result
         if (t >= static_cast<double>(nsave) * dt_save || t == Tf)
@@ -219,7 +218,7 @@ void run_simulation(Field& u,
             nsave++;
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 }
 
@@ -231,7 +230,7 @@ int main(int argc, char* argv[])
 
     auto& app = samurai::initialize("Finite volume example for the Burgers equation in 1d", argc, argv);
 
-    std::cout << "------------------------- Burgers -------------------------" << std::endl;
+    samurai::io::print("------------------------- Burgers -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -301,12 +300,12 @@ int main(int argc, char* argv[])
     scheme.enable_max_level_flux(true);
     run_simulation(u, unp1, u_max, unp1_max, scheme, cfl, mra_config, init_sol, nfiles, path, filename, nsave);
 
-    std::cout << std::endl;
-    std::cout << "Run the following commands to view the results:" << std::endl;
-    std::cout << "max-level-flux disabled:" << std::endl;
-    std::cout << "     python ../python/read_mesh.py " << filename << "_recons_ite_ --field u --start 0 --end " << nsave << std::endl;
-    std::cout << "max-level-flux enabled:" << std::endl;
-    std::cout << "     python ../python/read_mesh.py " << filename << "_mlf_recons_ite_ --field u --start 0 --end " << nsave << std::endl;
+    samurai::io::print("\n");
+    samurai::io::print("Run the following commands to view the results:\n");
+    samurai::io::print("max-level-flux disabled:\n");
+    samurai::io::print("     python ../python/read_mesh.py {}_recons_ite_ --field u --start 0 --end {}\n", filename, nsave);
+    samurai::io::print("max-level-flux enabled:\n");
+    samurai::io::print("     python ../python/read_mesh.py {}_mlf_recons_ite_ --field u --start 0 --end {}\n", filename, nsave);
 
     samurai::finalize();
     return 0;
diff --git a/demos/FiniteVolume/burgers_os.cpp b/demos/FiniteVolume/burgers_os.cpp
index e9965ec1c..553581af7 100644
--- a/demos/FiniteVolume/burgers_os.cpp
+++ b/demos/FiniteVolume/burgers_os.cpp
@@ -12,6 +12,7 @@
 namespace fs = std::filesystem;
 
 #include <numbers>
+#include <xtensor/xio.hpp>
 
 #include "convection_nonlin_os.hpp"
 
@@ -51,7 +52,7 @@ void check_diff(auto& field, auto& lca_left, auto& lca_right)
             set(
                 [&](auto& i, auto)
                 {
-                    std::cout << "Difference found !! " << level << " " << i << "\n";
+                    samurai::io::print("Difference found !! {} {}\n", level, i);
                     auto level_ = samurai::make_scalar_field<std::size_t>("level", mesh);
                     samurai::for_each_cell(mesh,
                                            [&](const auto& cell)
@@ -67,11 +68,11 @@ void check_diff(auto& field, auto& lca_left, auto& lca_right)
                 {
                     if (xt::any(xt::abs(field(level, i) - field(level, i + (1 << level))) > 1e-13))
                     {
-                        std::cout << fmt::format("\nDifference found at level {} on interval {}:\n", level, i);
-                        std::cout << fmt::format("\tleft = {}\n", field(level, i));
-                        std::cout << fmt::format("\tright = {}\n", field(level, i + (1 << level)));
-                        std::cout << fmt::format("\terror = {}\n", xt::abs(field(level, i) - field(level, i + (1 << level))));
-                        std::cout << mesh << std::endl;
+                        samurai::io::print("\nDifference found at level {} on interval {}:\n", level, i);
+                        samurai::io::print("\tleft = {}\n", fmt::streamed(field(level, i)));
+                        samurai::io::print("\tright = {}\n", fmt::streamed(field(level, i + (1 << level))));
+                        samurai::io::print("\terror = {}\n", fmt::streamed(xt::abs(field(level, i) - field(level, i + (1 << level)))));
+                        samurai::io::print("{}\n", fmt::streamed(mesh));
                         auto level_ = samurai::make_scalar_field<std::size_t>("level", mesh);
                         samurai::for_each_cell(mesh,
                                                [&](const auto& cell)
@@ -94,7 +95,7 @@ int main(int argc, char* argv[])
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
-    std::cout << "------------------------- Burgers 1D -------------------------" << std::endl;
+    samurai::io::print("------------------------- Burgers 1D -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -134,7 +135,7 @@ int main(int argc, char* argv[])
 
     SAMURAI_PARSE(argc, argv);
 
-    std::cout << "  max_level = " << max_level << "   min_level = " << min_level << std::endl;
+    samurai::io::print("  max_level = {}   min_level = {}\n", max_level, min_level);
 
     //--------------------//
     // Problem definition //
@@ -199,7 +200,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.12f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.12f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -210,7 +211,7 @@ int main(int argc, char* argv[])
         }
         catch (...)
         {
-            std::cout << "Exception caught in check_diff after adaptation" << std::endl;
+            samurai::io::print("Exception caught in check_diff after adaptation\n");
             samurai::finalize();
             return 1;
         }
@@ -228,7 +229,7 @@ int main(int argc, char* argv[])
         }
         catch (...)
         {
-            std::cout << "Exception caught in check_diff after integration" << std::endl;
+            samurai::io::print("Exception caught in check_diff after integration\n");
             samurai::finalize();
             return 1;
         }
@@ -236,13 +237,13 @@ int main(int argc, char* argv[])
         // Save the result
         if (nfiles == 0 || t >= static_cast<double>(nsave) * dt_save || t == Tf)
         {
-            std::cout << "  (saving results)" << std::flush;
+            samurai::io::print("  (saving results)");
             std::string suffix = (nfiles != 1) ? fmt::format("_level_{}_{}_ite_{}", min_level, max_level, nsave) : "";
             save(path, filename, u, suffix);
             nsave++;
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     samurai::finalize();
diff --git a/demos/FiniteVolume/heat.cpp b/demos/FiniteVolume/heat.cpp
index 0062ae1dc..5fc690966 100644
--- a/demos/FiniteVolume/heat.cpp
+++ b/demos/FiniteVolume/heat.cpp
@@ -5,6 +5,7 @@
 #include <samurai/io/restart.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -53,7 +54,7 @@ int main(int argc, char* argv[])
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
-    std::cout << "------------------------- Heat -------------------------" << std::endl;
+    samurai::io::print("------------------------- Heat -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -199,7 +200,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -232,18 +233,16 @@ int main(int argc, char* argv[])
                                              {
                                                  return exact_solution(coord, t, diff_coeff);
                                              });
-            std::cout.precision(2);
-            std::cout << ", L2-error: " << std::scientific << error;
+            samurai::io::print(", L2-error: {:.2e}", error);
         }
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     if (!save_final_state_only && dim == 1)
     {
-        std::cout << std::endl;
-        std::cout << "Run the following command to view the results:" << std::endl;
-        std::cout << "python <<path to samurai>>/python/read_mesh.py " << filename << "_ite_ --field u level --start 1 --end " << nsave
-                  << std::endl;
+        samurai::io::print("\n");
+        samurai::io::print("Run the following command to view the results:\n");
+        samurai::io::print("python <<path to samurai>>/python/read_mesh.py {}_ite_ --field u level --start 1 --end {}\n", filename, nsave);
     }
 
     if (save_final_state_only)
diff --git a/demos/FiniteVolume/heat_heterogeneous.cpp b/demos/FiniteVolume/heat_heterogeneous.cpp
index 8b4823da2..7269c3ce2 100644
--- a/demos/FiniteVolume/heat_heterogeneous.cpp
+++ b/demos/FiniteVolume/heat_heterogeneous.cpp
@@ -5,6 +5,7 @@
 #include <samurai/io/restart.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -41,7 +42,7 @@ int main(int argc, char* argv[])
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
-    std::cout << "------------------------- Heat -------------------------" << std::endl;
+    samurai::io::print("------------------------- Heat -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -176,7 +177,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -203,7 +204,7 @@ int main(int argc, char* argv[])
             save(path, filename, u, fmt::format("_ite_{}", nsave++));
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     if (save_final_state_only)
diff --git a/demos/FiniteVolume/heat_nonlinear.cpp b/demos/FiniteVolume/heat_nonlinear.cpp
index c9c259e11..c2e6d8e8e 100644
--- a/demos/FiniteVolume/heat_nonlinear.cpp
+++ b/demos/FiniteVolume/heat_nonlinear.cpp
@@ -5,6 +5,7 @@
 #include <samurai/io/restart.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/schemes/fv.hpp>
 
 #include <filesystem>
@@ -106,7 +107,7 @@ int main(int argc, char* argv[])
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
-    std::cout << "------------------------- Non-linear heat -------------------------" << std::endl;
+    samurai::io::print("------------------------- Non-linear heat -------------------------\n");
 
     /*
         Solves the non-linear heat equation
@@ -227,7 +228,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Update boundary conditions
         if (explicit_scheme)
@@ -271,8 +272,7 @@ int main(int argc, char* argv[])
                                          {
                                              return exact_solution(coords, t);
                                          });
-        std::cout.precision(2);
-        std::cout << ", L2-error: " << std::scientific << error;
+        samurai::io::print(", L2-error: {:.2e}", error);
 
         // Save the result
         if (!save_final_state_only)
@@ -280,7 +280,7 @@ int main(int argc, char* argv[])
             save(path, filename, u, fmt::format("_ite_{}", nsave++));
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     if (save_final_state_only)
diff --git a/demos/FiniteVolume/level_set.cpp b/demos/FiniteVolume/level_set.cpp
index 49d2b7c60..53198993e 100644
--- a/demos/FiniteVolume/level_set.cpp
+++ b/demos/FiniteVolume/level_set.cpp
@@ -331,7 +331,7 @@ int main(int argc, char* argv[])
         std::size_t ite = 0;
         while (true)
         {
-            std::cout << "Mesh adaptation iteration " << ite++ << std::endl;
+            samurai::io::print("Mesh adaptation iteration {}\n", ite++);
             tag.resize();
             AMR_criteria(phi, tag);
             samurai::graduation(tag, stencil_grad);
@@ -349,7 +349,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         // Numerical scheme
         samurai::update_ghost(phi, u);
diff --git a/demos/FiniteVolume/level_set_from_scratch.cpp b/demos/FiniteVolume/level_set_from_scratch.cpp
index 2000b1adc..19a614f99 100644
--- a/demos/FiniteVolume/level_set_from_scratch.cpp
+++ b/demos/FiniteVolume/level_set_from_scratch.cpp
@@ -667,7 +667,7 @@ int main(int argc, char* argv[])
         std::size_t ite = 0;
         while (true)
         {
-            std::cout << "Mesh adaptation iteration " << ite++ << std::endl;
+            samurai::io::print("Mesh adaptation iteration {}\n", ite++);
             auto tag = samurai::make_scalar_field<int>("tag", mesh);
             AMR_criteria(phi, tag);
             ::make_graduation(tag);
@@ -686,7 +686,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         // Numerical scheme
         update_ghosts(phi, u);
diff --git a/demos/FiniteVolume/lid_driven_cavity.cpp b/demos/FiniteVolume/lid_driven_cavity.cpp
index 004cdde63..201709b9c 100644
--- a/demos/FiniteVolume/lid_driven_cavity.cpp
+++ b/demos/FiniteVolume/lid_driven_cavity.cpp
@@ -148,7 +148,7 @@ int main(int argc, char* argv[])
 
     auto box = samurai::Box<double, dim>({0, 0}, {1, 1});
 
-    std::cout << "lid-driven cavity" << std::endl;
+    samurai::io::print("lid-driven cavity\n");
 
     //-------------------- 1 -----------------------------------------------------------------
     //
@@ -322,7 +322,7 @@ int main(int argc, char* argv[])
             t              = Tf;
             dt_has_changed = true;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt);
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation for Navier-Stokes
         if (min_level != max_level)
@@ -338,9 +338,9 @@ int main(int argc, char* argv[])
                 zero.resize();
                 rhs.resize();
             }
-            std::cout << ", levels " << min_level_np1 << "-" << max_level_np1;
+            samurai::io::print(", levels {}-{}", min_level_np1, max_level_np1);
         }
-        std::cout << std::endl;
+        samurai::io::print("\n");
 
         // Update solver
         if (mesh_has_changed || dt_has_changed)
@@ -443,20 +443,19 @@ int main(int argc, char* argv[])
         // Vec nest_result           = nestedAssembly.create_rhs_vector(result_nest);
         // MatMult(nestedA, nest_x, nest_result);
 
-        // std::cout << std::setprecision(15);
-        // std::cout << std::fixed;
+        // samurai::io::print("(setprecision 15, fixed)\n");
         // samurai::for_each_cell(
         //     mesh[mesh_id_t::reference],
         //     [&](auto cell)
         //     {
-        //         std::cout << round(result_velocity_mono[cell][0] * 1.e8) / 1.e8 << std::endl;
+        //         samurai::io::print("{}\n", round(result_velocity_mono[cell][0] * 1.e8) / 1.e8);
         //         if (round(result_velocity_mono[cell][0] * 1.e5) / 1.e5 != round(result_velocity_nest[cell][0] * 1.e5) / 1.e5)
         //         {
-        //             std::cout << result_velocity_mono[cell][0] << " =? " << result_velocity_nest[cell][0] << std::endl;
+        //             samurai::io::print("{} =? {}\n", result_velocity_mono[cell][0], result_velocity_nest[cell][0]);
         //         }
         //         if (result_pressure_mono[cell] != result_pressure_nest[cell])
         //         {
-        //             std::cout << result_pressure_mono[cell] << " =? " << result_pressure_nest[cell] << std::endl;
+        //             samurai::io::print("{} =? {}\n", result_pressure_mono[cell], result_pressure_nest[cell]);
         //         }
         //     });
     }
diff --git a/demos/FiniteVolume/linear_convection.cpp b/demos/FiniteVolume/linear_convection.cpp
index a7ae59070..cc432d67f 100644
--- a/demos/FiniteVolume/linear_convection.cpp
+++ b/demos/FiniteVolume/linear_convection.cpp
@@ -5,6 +5,7 @@
 #include <samurai/io/restart.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -41,7 +42,7 @@ int main(int argc, char* argv[])
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
-    std::cout << "------------------------- Linear convection -------------------------" << std::endl;
+    samurai::io::print("------------------------- Linear convection -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -171,7 +172,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -226,15 +227,14 @@ int main(int argc, char* argv[])
             }
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     if constexpr (dim == 1)
     {
-        std::cout << std::endl;
-        std::cout << "Run the following command to view the results:" << std::endl;
-        std::cout << "python <<path to samurai>>/python/read_mesh.py " << filename << "_ite_ --field u level --start 0 --end " << nsave
-                  << std::endl;
+        samurai::io::print("\n");
+        samurai::io::print("Run the following command to view the results:\n");
+        samurai::io::print("python <<path to samurai>>/python/read_mesh.py {}_ite_ --field u level --start 0 --end {}\n", filename, nsave);
     }
 
     samurai::finalize();
diff --git a/demos/FiniteVolume/linear_convection_obstacle.cpp b/demos/FiniteVolume/linear_convection_obstacle.cpp
index 00a048657..7588464d6 100644
--- a/demos/FiniteVolume/linear_convection_obstacle.cpp
+++ b/demos/FiniteVolume/linear_convection_obstacle.cpp
@@ -5,6 +5,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -39,7 +40,7 @@ int main(int argc, char* argv[])
     using Config                     = samurai::MRConfig<dim, 3>;
     using Mesh                       = samurai::MRMesh<Config>;
 
-    std::cout << "------------------------- Linear convection -------------------------" << std::endl;
+    samurai::io::print("------------------------- Linear convection -------------------------\n");
 
     //--------------------//
     // Program parameters //
@@ -145,7 +146,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
 
@@ -185,7 +186,7 @@ int main(int argc, char* argv[])
             }
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     samurai::finalize();
diff --git a/demos/FiniteVolume/manual_block_matrix_assembly.cpp b/demos/FiniteVolume/manual_block_matrix_assembly.cpp
index 03b86f243..18b7ec0be 100644
--- a/demos/FiniteVolume/manual_block_matrix_assembly.cpp
+++ b/demos/FiniteVolume/manual_block_matrix_assembly.cpp
@@ -4,6 +4,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -191,7 +192,7 @@ int main(int argc, char* argv[])
     using Config                     = samurai::MRConfig<dim>;
     using Box                        = samurai::Box<double, dim>;
 
-    std::cout << "------------------------- Begin -------------------------" << std::endl;
+    samurai::io::print("------------------------- Begin -------------------------\n");
 
     std::size_t min_level = 3;
     std::size_t max_level = 3;
@@ -275,18 +276,18 @@ int main(int argc, char* argv[])
     // Insert the coefficients into the matrix
     assembly.assemble_matrix(J);
 
-    std::cout << "Useless ghost rows: ";
+    samurai::io::print("Useless ghost rows: ");
     // assembly.get<0, 0>().for_each_useless_ghost_row(
     assembly.for_each_useless_ghost_row(
         [](auto row)
         {
-            std::cout << row << " ";
+            samurai::io::print("{} ", row);
         });
-    std::cout << std::endl;
+    samurai::io::print("\n");
 
     Vec v = assembly.create_vector(u_e, aux_Ce, u_s);
     VecView(v, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-    std::cout << std::endl;
+    samurai::io::print("\n");
 
     samurai::finalize();
     return 0;
diff --git a/demos/FiniteVolume/nagumo.cpp b/demos/FiniteVolume/nagumo.cpp
index 5f051d90e..80d25d314 100644
--- a/demos/FiniteVolume/nagumo.cpp
+++ b/demos/FiniteVolume/nagumo.cpp
@@ -4,6 +4,7 @@
 #include <samurai/io/restart.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -41,7 +42,7 @@ int main(int argc, char* argv[])
     using Box                           = samurai::Box<double, dim>;
     using point_t                       = typename Box::point_t;
 
-    std::cout << "------------------------- Nagumo -------------------------" << std::endl;
+    samurai::io::print("------------------------- Nagumo -------------------------\n");
 
     /**
      * Nagumo, or Fisher-KPP equation:
@@ -203,7 +204,7 @@ int main(int argc, char* argv[])
             dt += Tf - t;
             t = Tf;
         }
-        std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt) << std::flush;
+        samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt));
 
         // Mesh adaptation
         MRadaptation(mra_config);
@@ -252,8 +253,7 @@ int main(int argc, char* argv[])
                                          {
                                              return exact_solution(coord(0), t);
                                          });
-        std::cout.precision(2);
-        std::cout << ", L2-error: " << std::scientific << error;
+        samurai::io::print(", L2-error: {:.2e}", error);
 
         // Save the result
         if (!save_final_state_only)
@@ -261,15 +261,14 @@ int main(int argc, char* argv[])
             save(path, filename, u, fmt::format("_ite_{}", nsave++));
         }
 
-        std::cout << std::endl;
+        samurai::io::print("\n");
     }
 
     if (!save_final_state_only && dim == 1)
     {
-        std::cout << std::endl;
-        std::cout << "Run the following command to view the results:" << std::endl;
-        std::cout << "python <<path to samurai>>/python/read_mesh.py " << filename << "_ite_ --field u level --start 1 --end " << nsave
-                  << std::endl;
+        samurai::io::print("\n");
+        samurai::io::print("Run the following command to view the results:\n");
+        samurai::io::print("python <<path to samurai>>/python/read_mesh.py {}_ite_ --field u level --start 1 --end {}\n", filename, nsave);
     }
 
     if (save_final_state_only)
diff --git a/demos/FiniteVolume/scalar_burgers_2d.cpp b/demos/FiniteVolume/scalar_burgers_2d.cpp
index 0d062a215..11d21edf8 100644
--- a/demos/FiniteVolume/scalar_burgers_2d.cpp
+++ b/demos/FiniteVolume/scalar_burgers_2d.cpp
@@ -264,7 +264,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration {}: t = {}, dt = {}\n", nt++, t, dt);
 
         samurai::update_ghost_mr(u);
         unp1.resize();
diff --git a/demos/FiniteVolume/stencil_field.hpp b/demos/FiniteVolume/stencil_field.hpp
index 8adbf15ab..bd1841ed4 100644
--- a/demos/FiniteVolume/stencil_field.hpp
+++ b/demos/FiniteVolume/stencil_field.hpp
@@ -317,15 +317,15 @@ namespace samurai
         template <class T1>
         inline auto left_flux(const T1& u, double lb) const
         {
-            // std::cout << "left flux " << level << " " << i << " " << lb << std::endl;
-            // std::cout << flux(u(level, i - 1), u(level, i), lb) << std::endl;
+            // samurai::io::print("left flux {} {} {}\n", level, i, lb);
+            // samurai::io::print("{}\n", flux(u(level, i - 1), u(level, i), lb));
             return flux(u(level, i - 1), u(level, i), lb);
         }
 
         template <class T1>
         inline auto right_flux(const T1& u, double lb) const
         {
-            // std::cout << flux(u(level, i), u(level, i + 1), lb) << std::endl;
+            // samurai::io::print("{}\n", flux(u(level, i), u(level, i + 1), lb));
             return flux(u(level, i), u(level, i + 1), lb);
         }
     };
diff --git a/demos/FiniteVolume/stokes_2d.cpp b/demos/FiniteVolume/stokes_2d.cpp
index 463786eb7..68d2ba2bd 100644
--- a/demos/FiniteVolume/stokes_2d.cpp
+++ b/demos/FiniteVolume/stokes_2d.cpp
@@ -5,6 +5,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -21,7 +22,7 @@ template <class Field>
         if (std::isnan(value) || std::isinf(value) || (abs(value) < 1e-300 && abs(value) != 0))
         {
             is_nan_or_inf = true;
-            std::cout << f << std::endl;
+            samurai::io::print("{}\n", fmt::streamed(f));
             break;
         }
     }
@@ -203,11 +204,11 @@ int main(int argc, char* argv[])
     // Stationary problem //
     //--------------------//
 
-    std::cout << "Problem solved: ";
+    samurai::io::print("Problem solved: ");
 
     if (test_case == "s")
     {
-        std::cout << "stationary" << std::endl;
+        samurai::io::print("stationary\n");
         if (filename.empty())
         {
             filename = "stokes";
@@ -273,14 +274,14 @@ int main(int argc, char* argv[])
         zero.fill(0);
 
         // Linear solver
-        std::cout << "Solving Stokes system..." << std::endl;
+        samurai::io::print("Solving Stokes system...\n");
         auto stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);
 
         stokes_solver.set_unknowns(velocity, pressure);
         configure_solver(stokes_solver);
 
         stokes_solver.solve(f, zero);
-        std::cout << stokes_solver.iterations() << " iterations" << std::endl << std::endl;
+        samurai::io::print("{} iterations\n\n", stokes_solver.iterations());
 
         // Error
         double error = L2_error(velocity,
@@ -292,11 +293,10 @@ int main(int argc, char* argv[])
                                     auto v_y      = -v_x;
                                     return samurai::Array<double, dim, is_soa>{v_x, v_y};
                                 });
-        std::cout.precision(2);
-        std::cout << "L2-error on the velocity: " << std::scientific << error << std::endl;
+        samurai::io::print("L2-error on the velocity: {:.2e}\n", error);
 
         // Save solution
-        std::cout << "Saving solution..." << std::endl;
+        samurai::io::print("Saving solution...\n");
 
         samurai::save(path, filename, mesh, velocity);
         samurai::save(path, "pressure", mesh, pressure);
@@ -337,7 +337,7 @@ int main(int argc, char* argv[])
 
     else if (test_case == "ns")
     {
-        std::cout << "non stationary" << std::endl;
+        samurai::io::print("non stationary\n");
 
         if (filename.empty())
         {
@@ -464,7 +464,7 @@ int main(int argc, char* argv[])
                 t_np1          = Tf;
                 dt_has_changed = true;
             }
-            std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t_np1, dt);
+            samurai::io::print(samurai::io::root, "{}", fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t_np1, dt));
 
             // Mesh adaptation
             if (min_level != max_level)
@@ -480,9 +480,9 @@ int main(int argc, char* argv[])
                     rhs.resize();
                     zero.resize();
                 }
-                std::cout << ", levels " << min_level_np1 << "-" << max_level_np1;
+                samurai::io::print(", levels {}-{}", min_level_np1, max_level_np1);
             }
-            std::cout.flush();
+            // flush no-op with fmt
 
             // Boundary conditions
             velocity_np1.get_bc().clear();
@@ -536,8 +536,7 @@ int main(int argc, char* argv[])
                                     {
                                         return exact_velocity(t_n, coord);
                                     });
-            std::cout.precision(2);
-            std::cout << ", L2-error: " << std::scientific << error;
+            samurai::io::print(", L2-error: {:.2e}", error);
 
             // Divergence
             auto div_velocity = div(velocity);
@@ -558,8 +557,7 @@ int main(int argc, char* argv[])
                                                         {
                                                             return exact_velocity(t_n, coord);
                                                         });
-                std::cout.precision(2);
-                std::cout << ", L2-error (recons): " << std::scientific << error_recons;
+                samurai::io::print(", L2-error (recons): {:.2e}", error_recons);
                 // Save
                 if (nfiles != 1)
                 {
@@ -570,7 +568,7 @@ int main(int argc, char* argv[])
             {
                 nsave++;
             }
-            std::cout << std::endl;
+            samurai::io::print("\n");
         }
 
         if (nfiles == 1)
@@ -580,7 +578,7 @@ int main(int argc, char* argv[])
     }
     else
     {
-        std::cerr << "Unknown test case. Allowed options are 's' = stationary, 'ns' = non-stationary." << std::endl;
+        samurai::io::eprint("Unknown test case. Allowed options are 's' = stationary, 'ns' = non-stationary.\n");
     }
 
     samurai::finalize();
diff --git a/demos/LBM/D1Q222_Euler_Sod.cpp b/demos/LBM/D1Q222_Euler_Sod.cpp
index a8622eb20..3c53cf86a 100644
--- a/demos/LBM/D1Q222_Euler_Sod.cpp
+++ b/demos/LBM/D1Q222_Euler_Sod.cpp
@@ -13,6 +13,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/subset/node.hpp>
 
@@ -470,7 +471,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -514,21 +515,20 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite < N; ++nb_ite)
             {
-                std::cout << std::endl << "Iteration = " << nb_ite << std::endl;
+                samurai::io::print("\nIteration = {}\n", nb_ite);
                 MRadaptation(eps, regularity);
 
                 save_solution(f, eps, nb_ite);
 
                 auto error = compute_error(f, fR, update_bc_for_level, t);
 
-                std::cout << std::endl
-                          << "Error rho = " << error[0] << std::endl
-                          << "Diff rho = " << error[1] << std::endl
-                          << "Error q = " << error[2] << std::endl
-                          << "Diff q = " << error[3] << std::endl
-                          << "Error E = " << error[4] << std::endl
-                          << "Diff E = " << error[5];
-                std::cout << std::endl;
+                samurai::io::print("\nError rho = {}\nDiff rho = {}\nError q = {}\nDiff q = {}\nError E = {}\nDiff E = {}\n",
+                                   error[0],
+                                   error[1],
+                                   error[2],
+                                   error[3],
+                                   error[4],
+                                   error[5]);
 
                 one_time_step(f, pred_coeff_separate, update_bc_for_level, s, lambda);
                 one_time_step(fR, pred_coeff_separate, update_bc_for_level, s, lambda);
@@ -538,7 +538,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/D1Q2_Advection_and_Burgers.cpp b/demos/LBM/D1Q2_Advection_and_Burgers.cpp
index 153b91fdd..a7a56af94 100644
--- a/demos/LBM/D1Q2_Advection_and_Burgers.cpp
+++ b/demos/LBM/D1Q2_Advection_and_Burgers.cpp
@@ -13,6 +13,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -134,7 +135,7 @@ xt::xtensor<double, 1> prediction(const Field& f,
     auto it = mem_map.find({item, level_g, level, i});
     if (it != mem_map.end())
     {
-        // std::cout<<std::endl<<"Found by memoization";
+        // samurai::io::print("\nFound by memoization");
         return it->second;
     }
     else
@@ -144,7 +145,7 @@ xt::xtensor<double, 1> prediction(const Field& f,
         xt::xtensor<double, 1> out = xt::empty<double>({i.size() / i.step}); // xt::eval(f(item, level_g, i));
         auto mask                  = mesh.exists(mesh_id_t::cells_and_ghosts, level_g + level, i);
 
-        // std::cout << level_g + level << " " << i << " " << mask << "\n";
+        // samurai::io::print("{} {} {}\n", level_g + level, i, mask);
         if (xt::all(mask))
         {
             return xt::eval(f(item, level_g + level, i));
@@ -534,7 +535,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -594,8 +595,7 @@ int main(int argc, char* argv[])
                 save_reconstructed(f, fR, update_bc_for_level, eps, nb_ite);
 
                 auto error = compute_error(f, fR, update_bc_for_level, t);
-
-                std::cout << std::endl << "Diff = " << error[1] << std::flush;
+                samurai::io::print("\nDiff = {}", error[1]);
 
                 out_time_frames << t << std::endl;
                 out_error_exact_ref << error[0] << std::endl;
@@ -618,7 +618,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/D1Q3_Shallow_Waters_Dam.cpp b/demos/LBM/D1Q3_Shallow_Waters_Dam.cpp
index 43207281e..030771dfc 100644
--- a/demos/LBM/D1Q3_Shallow_Waters_Dam.cpp
+++ b/demos/LBM/D1Q3_Shallow_Waters_Dam.cpp
@@ -13,6 +13,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -346,7 +347,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -393,18 +394,14 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite < N; ++nb_ite)
             {
-                std::cout << std::endl << "Iteration " << nb_ite << " Time = " << t;
+                samurai::io::print("\nIteration {} Time = {}", nb_ite, t);
 
                 MRadaptation(eps, 0.); // Regularity 0
 
                 save_solution(f, eps, nb_ite, lambda);
 
                 auto error = compute_error(f, fR, update_bc_for_level, t, lambda, g);
-                std::cout << std::endl
-                          << "Error h = " << error[0] << std::endl
-                          << "Diff h = " << error[1] << std::endl
-                          << "Error q = " << error[2] << std::endl
-                          << "Diff q = " << error[3];
+                samurai::io::print("\nError h = {}\nDiff h = {}\nError q = {}\nDiff q = {}", error[0], error[1], error[2], error[3]);
 
                 tic();
                 one_time_step(f, pred_coeff_separate, update_bc_for_level, s, lambda, g);
@@ -417,9 +414,9 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
-    std::cout << std::endl;
+    samurai::io::print("\n");
     samurai::finalize();
     return 0;
 }
diff --git a/demos/LBM/D1Q5_Shallow_Waters_Dam.cpp b/demos/LBM/D1Q5_Shallow_Waters_Dam.cpp
index eb51fe094..42b682beb 100644
--- a/demos/LBM/D1Q5_Shallow_Waters_Dam.cpp
+++ b/demos/LBM/D1Q5_Shallow_Waters_Dam.cpp
@@ -13,6 +13,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -583,7 +584,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -629,19 +630,14 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite < N; ++nb_ite)
             {
-                std::cout << std::endl << "Iteration " << nb_ite << " Time = " << t;
+                samurai::io::print("\nIteration {} Time = {}", nb_ite, t);
 
                 MRadaptation(eps, 0.); // Regularity 0
 
                 save_solution(f, eps, nb_ite, lambda);
 
                 auto error = compute_error(f, fR, update_bc_for_level, t, lambda, g);
-
-                std::cout << std::endl
-                          << "Error h = " << error[0] << std::endl
-                          << "Diff h = " << error[1] << std::endl
-                          << "Error q = " << error[2] << std::endl
-                          << "Diff q = " << error[3];
+                samurai::io::print("\nError h = {}\nDiff h = {}\nError q = {}\nDiff q = {}", error[0], error[1], error[2], error[3]);
 
                 tic();
                 one_time_step_overleaves(f, pred_coeff_separate, update_bc_for_level, s, lambda, g);
@@ -654,9 +650,9 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
-    std::cout << std::endl;
+    samurai::io::print("\n");
     samurai::finalize();
     return 0;
 }
diff --git a/demos/LBM/D2Q4444_Euler_Implosion.cpp b/demos/LBM/D2Q4444_Euler_Implosion.cpp
index 598f0ba12..8efcadd3e 100644
--- a/demos/LBM/D2Q4444_Euler_Implosion.cpp
+++ b/demos/LBM/D2Q4444_Euler_Implosion.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -391,7 +392,7 @@ void one_time_step(Field& f,
             std::size_t j = max_level - (lev_p_1);
             double coeff  = 1. / (1 << (2 * j)); // ATTENTION A LA DIMENSION 2 !!!!
 
-            // std::cout<<std::endl<<"Level = "<<level<<std::endl;
+            // samurai::io::print("\nLevel = {}\n", level);
 
             auto leaves = samurai::intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::cells][level]);
 
@@ -446,7 +447,7 @@ void one_time_step(Field& f,
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
 
-                    // std::cout<<std::endl<<"East : "<<k<<"  | 
+                    // samurai::io::print("\nEast : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{0, 2, 6};
@@ -479,7 +480,7 @@ void one_time_step(Field& f,
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
 
-                    // std::cout<<std::endl<<"NorthEast : "<<k<<"  | 
+                    // samurai::io::print("\nNorthEast : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{0, 2};
@@ -516,7 +517,7 @@ void one_time_step(Field& f,
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
 
-                    // std::cout<<std::endl<<"SouthEast : "<<k<<"  | 
+                    // samurai::io::print("\nSouthEast : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{0, 6};
@@ -552,7 +553,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"West : "<<k<<"  | 
+                    // samurai::io::print("\nWest : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{2, 4, 6};
@@ -584,7 +585,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"NorthWest : "<<k<<"  | 
+                    // samurai::io::print("\nNorthWest : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{2, 4};
@@ -620,7 +621,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"SouthWest : "<<k<<"  | 
+                    // samurai::io::print("\nSouthWest : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{4, 6};
@@ -656,7 +657,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"North : "<<k<<"  | 
+                    // samurai::io::print("\nNorth : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{0, 2, 4};
@@ -688,7 +689,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"South : "<<k<<"  | 
+                    // samurai::io::print("\nSouth : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{0, 4, 6};
@@ -1024,7 +1025,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -1093,7 +1094,7 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite <= N; ++nb_ite)
             {
-                std::cout << std::endl << "Time = " << t << "   Iteration number = " << nb_ite << std::endl;
+                samurai::io::print("\nTime = {}   Iteration number = {}\n", t, nb_ite);
 
                 if (max_level > min_level)
                 {
@@ -1114,7 +1115,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/D2Q4444_Euler_Lax_Liu.cpp b/demos/LBM/D2Q4444_Euler_Lax_Liu.cpp
index d45c07f86..7b050c791 100644
--- a/demos/LBM/D2Q4444_Euler_Lax_Liu.cpp
+++ b/demos/LBM/D2Q4444_Euler_Lax_Liu.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/statistics.hpp>
 
@@ -779,7 +780,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -880,7 +881,7 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite <= N; ++nb_ite)
             {
-                std::cout << std::endl << "   Iteration number = " << nb_ite << std::endl;
+                samurai::io::print("\n   Iteration number = {}\n", nb_ite);
 
                 if (max_level > min_level)
                 {
@@ -890,7 +891,7 @@ int main(int argc, char* argv[])
                 if (nb_ite == N)
                 {
                     auto error_density = compute_error(f, f_ref, update_bc_for_level);
-                    std::cout << std::endl << "#### Epsilon = " << eps << "   error = " << error_density << std::endl;
+                    samurai::io::print("\n#### Epsilon = {}   error = {}\n", eps, error_density);
                     save_solution(f, eps, nb_ite, std::string("final_"));
                     save_reconstructed(f, f_ref, update_bc_for_level, eps, nb_ite);
                 }
@@ -923,7 +924,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/D2Q4444_Euler_Lax_Liu_uniform.cpp b/demos/LBM/D2Q4444_Euler_Lax_Liu_uniform.cpp
index d58a65e49..70653da0d 100644
--- a/demos/LBM/D2Q4444_Euler_Lax_Liu_uniform.cpp
+++ b/demos/LBM/D2Q4444_Euler_Lax_Liu_uniform.cpp
@@ -8,6 +8,7 @@
 
 #include <samurai/field.hpp>
 #include <samurai/io/hdf5.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/statistics.hpp>
 #include <samurai/uniform_mesh.hpp>
@@ -355,7 +356,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -414,7 +415,7 @@ int main(int argc, char* argv[])
             tic();
             for (std::size_t nb_ite = 0; nb_ite <= N; ++nb_ite)
             {
-                // std::cout<<std::endl<<"   Iteration number =
+                // samurai::io::print("\n   Iteration number = ");
                 // "<<nb_ite<<std::endl;
 
                 if (nb_ite == N)
@@ -427,13 +428,12 @@ int main(int argc, char* argv[])
                 // samurai::statistics("D2Q4444_Euler_Lax_Liu", mesh);
             }
             auto duration = toc();
-            std::cout << "nb_ite: " << N << " execution time: " << duration
-                      << " MLUPS: " << N * (1 << level) * (1 << level) / duration / 1e6 << std::endl;
+            samurai::io::print("nb_ite: {} execution time: {} MLUPS: {}\n", N, duration, N * (1 << level) * (1 << level) / duration / 1e6);
         }
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/D2Q5444_Euler_Rayleigh_Taylor.cpp b/demos/LBM/D2Q5444_Euler_Rayleigh_Taylor.cpp
index de8e8b846..2cd65eb9e 100644
--- a/demos/LBM/D2Q5444_Euler_Rayleigh_Taylor.cpp
+++ b/demos/LBM/D2Q5444_Euler_Rayleigh_Taylor.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -510,7 +511,7 @@ void one_time_step(Field& f,
             std::size_t j = max_level - (lev_p_1);
             double coeff  = 1. / (1 << (2 * j)); // ATTENTION A LA DIMENSION 2 !!!!
 
-            // std::cout<<std::endl<<"Level = "<<level<<std::endl;
+            // samurai::io::print("\nLevel = {}\n", level);
 
             auto leaves = samurai::intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::cells][level]);
 
@@ -565,7 +566,7 @@ void one_time_step(Field& f,
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
 
-                    // std::cout<<std::endl<<"East : "<<k<<"  | 
+                    // samurai::io::print("\nEast : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{0, 2, 6};
@@ -604,7 +605,7 @@ void one_time_step(Field& f,
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
 
-                    // std::cout<<std::endl<<"NorthEast : "<<k<<"  | 
+                    // samurai::io::print("\nNorthEast : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{0, 2};
@@ -648,7 +649,7 @@ void one_time_step(Field& f,
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
 
-                    // std::cout<<std::endl<<"SouthEast : "<<k<<"  | 
+                    // samurai::io::print("\nSouthEast : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{0, 6};
@@ -689,7 +690,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"West : "<<k<<"  | 
+                    // samurai::io::print("\nWest : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{2, 4, 6};
@@ -727,7 +728,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"NorthWest : "<<k<<"  | 
+                    // samurai::io::print("\nNorthWest : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{2, 4};
@@ -769,7 +770,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"SouthWest : "<<k<<"  | 
+                    // samurai::io::print("\nSouthWest : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 2> vld_flx{4, 6};
@@ -810,7 +811,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"North : "<<k<<"  | 
+                    // samurai::io::print("\nNorth : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{0, 2, 4};
@@ -849,7 +850,7 @@ void one_time_step(Field& f,
                 {
                     auto k = interval; // Logical index in x
                     auto h = index[0]; // Logical index in y
-                    // std::cout<<std::endl<<"South : "<<k<<"  | 
+                    // samurai::io::print("\nSouth : {}  | ", k);
                     // "<<h<<std::endl;
 
                     std::array<unsigned short int, 3> vld_flx{0, 4, 6};
@@ -1104,7 +1105,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -1152,7 +1153,7 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite <= N; ++nb_ite)
             {
-                std::cout << std::endl << "Time = " << t << "   Iteration number = " << nb_ite << std::endl;
+                samurai::io::print("\nTime = {}   Iteration number = {}\n", t, nb_ite);
 
                 if (max_level > min_level)
                 {
@@ -1173,7 +1174,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/D2Q9_Navier_Stokes_von_Karman_street.cpp b/demos/LBM/D2Q9_Navier_Stokes_von_Karman_street.cpp
index 26df36daa..88d30afa4 100644
--- a/demos/LBM/D2Q9_Navier_Stokes_von_Karman_street.cpp
+++ b/demos/LBM/D2Q9_Navier_Stokes_von_Karman_street.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -1359,7 +1360,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print("{}\n", options.help());
         }
         else
         {
@@ -1424,7 +1425,7 @@ int main(int argc, char* argv[])
 
             for (std::size_t nb_ite = 0; nb_ite < N; ++nb_ite)
             {
-                std::cout << std::endl << "Iteration number = " << nb_ite << std::endl;
+                samurai::io::print("\nIteration number = {}\n", nb_ite);
                 time_frames << nb_ite * dt << std::endl;
 
                 if (max_level > min_level)
@@ -1440,7 +1441,7 @@ int main(int argc, char* argv[])
                 }
 
                 auto CDCL = one_time_step(f, update_bc_for_level, pred_coeff, rho0, u0, lambda, mu, zeta, radius, momenti);
-                std::cout << std::endl << "CD = " << CDCL.first << "   CL = " << CDCL.second << std::endl;
+                samurai::io::print("\nCD = {}   CL = {}\n", CDCL.first, CDCL.second);
                 CD << CDCL.first << std::endl;
                 CL << CDCL.second << std::endl;
 
@@ -1457,7 +1458,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print("{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/prediction_map_1d.hpp b/demos/LBM/prediction_map_1d.hpp
index 5aee55928..b16f3bc55 100644
--- a/demos/LBM/prediction_map_1d.hpp
+++ b/demos/LBM/prediction_map_1d.hpp
@@ -72,7 +72,7 @@ class prediction_map
     {
         for (auto& c : coeff)
         {
-            std::cout << "( " << c.first << ", ): " << c.second << "\n";
+            out << "( " << c.first << ", ): " << c.second << "\n";
         }
     }
 
diff --git a/demos/LBM/prediction_map_2d.hpp b/demos/LBM/prediction_map_2d.hpp
index 1c5354532..2e47b1213 100644
--- a/demos/LBM/prediction_map_2d.hpp
+++ b/demos/LBM/prediction_map_2d.hpp
@@ -72,7 +72,7 @@ class prediction_map
     {
         for (auto& c : coeff)
         {
-            std::cout << "( " << std::get<0>(c.first) << ", " << std::get<1>(c.first) << " ): " << c.second << "\n";
+            out << "( " << std::get<0>(c.first) << ", " << std::get<1>(c.first) << " ): " << c.second << "\n";
         }
     }
 
diff --git a/demos/LBM/test_D1Q2.cpp b/demos/LBM/test_D1Q2.cpp
index 1ea703f47..090d73be0 100644
--- a/demos/LBM/test_D1Q2.cpp
+++ b/demos/LBM/test_D1Q2.cpp
@@ -11,6 +11,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/reconstruction.hpp>
 #include <samurai/samurai.hpp>
 
@@ -369,7 +370,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print(samurai::io::root, "{}\n", options.help());
         }
         else
         {
@@ -445,11 +446,11 @@ int main(int argc, char* argv[])
 
             for (auto s : s_vect)
             {
-                std::cout << std::endl << "Relaxation parameter s = " << s;
+                samurai::io::print(samurai::io::root, "\nRelaxation parameter s = {}", s);
 
                 std::string prefix(case_name + "_s_" + std::to_string(s) + "_");
 
-                std::cout << std::endl << "Testing time behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting time behavior\n");
                 {
                     double eps = 1.0e-4; // This remains fixed
 
@@ -492,7 +493,7 @@ int main(int argc, char* argv[])
                                                / static_cast<double>(meshR.nb_cells(mesh_id_t::cells))
                                         << std::endl;
 
-                        std::cout << std::endl << "n = " << nb_ite << "   Time = " << t << " Diff = " << error[1] << std::endl;
+                        samurai::io::print(samurai::io::root, "\nn = {}   Time = {} Diff = {}\n", nb_ite, t, error[1]);
                         ;
 
                         one_time_step(f, update_bc_for_level, s, lambda, ad_vel, test_number, finest_collision);
@@ -500,7 +501,7 @@ int main(int argc, char* argv[])
                         t += dt;
                     }
 
-                    std::cout << std::endl;
+                    samurai::io::print(samurai::io::root, "\n");
 
                     out_time_frames.close();
                     out_error_exact_ref.close();
@@ -508,7 +509,7 @@ int main(int argc, char* argv[])
                     out_compression.close();
                 }
 
-                std::cout << std::endl << "Testing eps behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting eps behavior\n");
                 {
                     double eps         = 0.1;
                     std::size_t N_test = 50;
@@ -525,7 +526,7 @@ int main(int argc, char* argv[])
 
                     for (std::size_t n_test = 0; n_test < N_test; ++n_test)
                     {
-                        std::cout << std::endl << "Test " << n_test << " eps = " << eps;
+                        samurai::io::print(samurai::io::root, "\nTest {} eps = {}", n_test, eps);
 
                         mesh_t mesh{box, min_level, max_level};
                         mesh_t meshR{box, max_level, max_level}; // This is the reference scheme
@@ -555,7 +556,7 @@ int main(int argc, char* argv[])
                         }
 
                         auto error = compute_error(f, fR, update_bc_for_level, t, ad_vel, test_number);
-                        std::cout << "Diff = " << error[1] << std::endl;
+                        samurai::io::print(samurai::io::root, "Diff = {}\n", error[1]);
 
                         std::size_t max_level_effective = mesh.min_level();
 
@@ -589,7 +590,7 @@ int main(int argc, char* argv[])
 
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print(samurai::io::root, "{}\n", options.help());
     }
 
     samurai::finalize();
diff --git a/demos/LBM/test_D1Q222.cpp b/demos/LBM/test_D1Q222.cpp
index e1ee32b20..3241c4e72 100644
--- a/demos/LBM/test_D1Q222.cpp
+++ b/demos/LBM/test_D1Q222.cpp
@@ -13,6 +13,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -440,7 +441,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print(samurai::io::root, "{}\n", options.help());
         }
         else
         {
@@ -474,11 +475,11 @@ int main(int argc, char* argv[])
 
             for (auto s : s_vect)
             {
-                std::cout << std::endl << "Relaxation parameter s = " << s;
+                samurai::io::print(samurai::io::root, "\nRelaxation parameter s = {}", s);
 
                 std::string prefix(case_name + "_s_" + std::to_string(s) + "_");
 
-                std::cout << std::endl << "Testing time behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting time behavior\n");
                 {
                     double eps = 1.0e-4; // This remains fixed
 
@@ -535,17 +536,17 @@ int main(int argc, char* argv[])
                                                / static_cast<double>(meshR.nb_cells(mesh_id_t::cells))
                                         << std::endl;
 
-                        std::cout << std::endl
-                                  << "Time = " << t << " Diff_h = " << error[1] << std::endl
-                                  << "Diff q = " << error[3] << std::endl
-                                  << "Diff E = " << error[5];
+                        samurai::io::print(samurai::io::root, "\n");
+                        << "Time = " << t << " Diff_h = " << error[1] << std::endl
+                        << "Diff q = " << error[3] << std::endl
+                        << "Diff E = " << error[5];
 
                         one_time_step(f, pred_coeff_separate, update_bc_for_level, s, lambda);
                         one_time_step(fR, pred_coeff_separate, update_bc_for_level, s, lambda);
                         t += dt;
                     }
 
-                    std::cout << std::endl;
+                    samurai::io::print(samurai::io::root, "\n");
 
                     out_time_frames.close();
                     out_error_rho_exact_ref.close();
@@ -557,7 +558,7 @@ int main(int argc, char* argv[])
                     out_compression.close();
                 }
 
-                std::cout << std::endl << "Testing eps behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting eps behavior\n");
                 {
                     double eps         = 1.0e-1; // 0.1;
                     std::size_t N_test = 50;     // 50;
@@ -577,7 +578,7 @@ int main(int argc, char* argv[])
 
                     for (std::size_t n_test = 0; n_test < N_test; ++n_test)
                     {
-                        std::cout << std::endl << "Test " << n_test << " eps = " << eps;
+                        samurai::io::print(samurai::io::root, "\nTest {} eps = {}", n_test, eps);
 
                         samurai::MROMesh<Config> mesh{box, min_level, max_level};
                         samurai::MROMesh<Config> meshR{box, max_level, max_level}; // This is the reference scheme
@@ -604,9 +605,8 @@ int main(int argc, char* argv[])
                         }
 
                         auto error = compute_error(f, fR, update_bc_for_level, t);
-                        std::cout << "Diff  h= " << error[1] << std::endl
-                                  << "Diff q = " << error[3] << std::endl
-                                  << "Diff E = " << error[5] << std::endl;
+                        samurai::io::print(samurai::io::root, "Diff  h= {}\n", error[1]);
+                        << "Diff q = " << error[3] << std::endl << "Diff E = " << error[5] << std::endl;
 
                         out_eps << eps << std::endl;
                         out_diff_rho_ref_adap << error[1] << std::endl;
@@ -631,7 +631,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print(samurai::io::root, "{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/LBM/test_D1Q3.cpp b/demos/LBM/test_D1Q3.cpp
index 970b65d3a..8dd553997 100644
--- a/demos/LBM/test_D1Q3.cpp
+++ b/demos/LBM/test_D1Q3.cpp
@@ -12,6 +12,7 @@
 #include <samurai/field.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -318,7 +319,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print(samurai::io::root, "{}\n", options.help());
         }
         else
         {
@@ -354,11 +355,11 @@ int main(int argc, char* argv[])
 
             for (auto s : s_vect)
             {
-                std::cout << std::endl << "Relaxation parameter s = " << s;
+                samurai::io::print(samurai::io::root, "\nRelaxation parameter s = {}", s);
 
                 std::string prefix(case_name + "_s_" + std::to_string(s) + "_");
 
-                std::cout << std::endl << "Testing time behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting time behavior\n");
                 {
                     double eps = 1.0e-4; // This remains fixed
 
@@ -409,14 +410,14 @@ int main(int argc, char* argv[])
                                                / static_cast<double>(meshR.nb_cells(mesh_id_t::cells))
                                         << std::endl;
 
-                        std::cout << std::endl << "Time = " << t << " Diff_h = " << error[1] << std::endl << "Diff q = " << error[3];
+                        samurai::io::print(samurai::io::root, "\nTime = {} Diff_h = {}\nDiff q = {}", t, error[1], error[3]);
 
                         one_time_step(f, pred_coeff_separate, update_bc_for_level, s, lambda, g);
                         one_time_step(fR, pred_coeff_separate, update_bc_for_level, s, lambda, g);
                         t += dt;
                     }
 
-                    std::cout << std::endl;
+                    samurai::io::print(samurai::io::root, "\n");
 
                     out_time_frames.close();
                     out_error_h_exact_ref.close();
@@ -426,7 +427,7 @@ int main(int argc, char* argv[])
                     out_compression.close();
                 }
 
-                std::cout << std::endl << "Testing eps behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting eps behavior\n");
                 {
                     double eps         = 1.0e-1; // 0.1;
                     std::size_t N_test = 50;     // 50;
@@ -443,7 +444,7 @@ int main(int argc, char* argv[])
 
                     for (std::size_t n_test = 0; n_test < N_test; ++n_test)
                     {
-                        std::cout << std::endl << "Test " << n_test << " eps = " << eps;
+                        samurai::io::print(samurai::io::root, "\nTest {} eps = {}", n_test, eps);
 
                         samurai::MROMesh<Config> mesh{box, min_level, max_level};
                         samurai::MROMesh<Config> meshR{box, max_level, max_level}; // This is the reference scheme
@@ -470,7 +471,7 @@ int main(int argc, char* argv[])
                         }
 
                         auto error = compute_error(f, fR, update_bc_for_level, t, lambda, g);
-                        std::cout << "Diff  h= " << error[1] << std::endl << "Diff q = " << error[3] << std::endl;
+                        samurai::io::print(samurai::io::root, "Diff  h= {}\nDiff q = {}\n", error[1], error[3]);
 
                         out_eps << eps << std::endl;
                         out_diff_h_ref_adap << error[1] << std::endl;
@@ -493,7 +494,7 @@ int main(int argc, char* argv[])
 
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print(samurai::io::root, "{}\n", options.help());
     }
 
     samurai::finalize();
diff --git a/demos/LBM/test_D1Q5.cpp b/demos/LBM/test_D1Q5.cpp
index 7b1c2682e..6a00b3ee4 100644
--- a/demos/LBM/test_D1Q5.cpp
+++ b/demos/LBM/test_D1Q5.cpp
@@ -13,6 +13,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh_with_overleaves.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "boundary_conditions.hpp"
@@ -554,7 +555,7 @@ int main(int argc, char* argv[])
 
         if (result.count("help"))
         {
-            std::cout << options.help() << "\n";
+            samurai::io::print(samurai::io::root, "{}\n", options.help());
         }
         else
         {
@@ -591,11 +592,11 @@ int main(int argc, char* argv[])
 
             for (auto s : s_vect)
             {
-                std::cout << std::endl << "Relaxation parameter s = " << s;
+                samurai::io::print(samurai::io::root, "\nRelaxation parameter s = {}", s);
 
                 std::string prefix(case_name + "_s_" + std::to_string(s) + "_");
 
-                std::cout << std::endl << "Testing time behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting time behavior\n");
                 {
                     double eps = 1.0e-4; // This remains fixed
 
@@ -646,14 +647,14 @@ int main(int argc, char* argv[])
                                                / static_cast<double>(meshR.nb_cells(mesh_id_t::cells))
                                         << std::endl;
 
-                        std::cout << std::endl << "Time = " << t << " Diff_h = " << error[1] << std::endl << "Diff q = " << error[3];
+                        samurai::io::print(samurai::io::root, "\nTime = {} Diff_h = {}\nDiff q = {}", t, error[1], error[3]);
 
                         one_time_step_overleaves(f, pred_coeff_separate, update_bc_for_level, s, lambda, g);
                         one_time_step_overleaves(fR, pred_coeff_separate, update_bc_for_level, s, lambda, g);
                         t += dt;
                     }
 
-                    std::cout << std::endl;
+                    samurai::io::print(samurai::io::root, "\n");
 
                     out_time_frames.close();
                     out_error_h_exact_ref.close();
@@ -663,7 +664,7 @@ int main(int argc, char* argv[])
                     out_compression.close();
                 }
 
-                std::cout << std::endl << "Testing eps behavior" << std::endl;
+                samurai::io::print(samurai::io::root, "\nTesting eps behavior\n");
                 {
                     double eps         = 1.0e-1; // 0.1;
                     std::size_t N_test = 50;     // 50;
@@ -680,7 +681,7 @@ int main(int argc, char* argv[])
 
                     for (std::size_t n_test = 0; n_test < N_test; ++n_test)
                     {
-                        std::cout << std::endl << "Test " << n_test << " eps = " << eps;
+                        samurai::io::print(samurai::io::root, "\nTest {} eps = {}", n_test, eps);
 
                         samurai::MROMesh<Config> mesh{box, min_level, max_level};
                         samurai::MROMesh<Config> meshR{box, max_level, max_level}; // This is the reference scheme
@@ -707,7 +708,7 @@ int main(int argc, char* argv[])
                         }
 
                         auto error = compute_error(f, fR, update_bc_for_level, t, lambda, g);
-                        std::cout << "Diff  h= " << error[1] << std::endl << "Diff q = " << error[3] << std::endl;
+                        samurai::io::print(samurai::io::root, "Diff  h= {}\nDiff q = {}\n", error[1], error[3]);
 
                         out_eps << eps << std::endl;
                         out_diff_h_ref_adap << error[1] << std::endl;
@@ -729,7 +730,7 @@ int main(int argc, char* argv[])
     }
     catch (const cxxopts::OptionException& e)
     {
-        std::cout << options.help() << "\n";
+        samurai::io::print(samurai::io::root, "{}\n", options.help());
     }
     samurai::finalize();
     return 0;
diff --git a/demos/Weno/VF_level_set_houc5_amr.cpp b/demos/Weno/VF_level_set_houc5_amr.cpp
index b6ae884a7..f9f3fdc6c 100644
--- a/demos/Weno/VF_level_set_houc5_amr.cpp
+++ b/demos/Weno/VF_level_set_houc5_amr.cpp
@@ -1,6 +1,7 @@
 // #include <spdlog/spdlog.h>
 #include <fmt/format.h>
 #include <fmt/ostream.h>
+#include <samurai/print.hpp>
 #include <xtensor/xfixed.hpp>
 #include <xtensor/xmasked_view.hpp>
 
@@ -290,12 +291,12 @@ int main()
     std::array<double, dim> a{1, 0};
     for (std::size_t ite = 0; ite < max_ite; ++ite)
     {
-        std::cout << "iteration: " << ite << std::endl;
+        samurai::io::print(samurai::io::root, "iteration: {}\n", ite);
 
         std::size_t ite_adapt = 0;
         while (1)
         {
-            std::cout << "\tmesh adaptation: " << ite_adapt++ << std::endl;
+            samurai::io::print(samurai::io::root, "\tmesh adaptation: {}\n", ite_adapt++);
             samurai::update_ghost(update_bc, field);
             auto tag = samurai::make_scalar_field<int>("tag", mesh);
             AMR_criteria(field, tag);
diff --git a/demos/Weno/VF_level_set_os5_amr.cpp b/demos/Weno/VF_level_set_os5_amr.cpp
index 41d1f4a56..d36888814 100644
--- a/demos/Weno/VF_level_set_os5_amr.cpp
+++ b/demos/Weno/VF_level_set_os5_amr.cpp
@@ -1,6 +1,7 @@
 // #include <spdlog/spdlog.h>
 #include <fmt/format.h>
 #include <fmt/ostream.h>
+#include <samurai/print.hpp>
 #include <xtensor/xfixed.hpp>
 #include <xtensor/xmasked_view.hpp>
 
@@ -309,12 +310,12 @@ int main()
     std::array<double, dim> a{1, 0};
     for (std::size_t ite = 0; ite < max_ite; ++ite)
     {
-        std::cout << "iteration: " << ite << std::endl;
+        samurai::io::print(samurai::io::root, "iteration: {}\n", ite);
 
         std::size_t ite_adapt = 0;
         while (1)
         {
-            std::cout << "\tmesh adaptation: " << ite_adapt++ << std::endl;
+            samurai::io::print(samurai::io::root, "\tmesh adaptation: {}\n", ite_adapt++);
             samurai::update_ghost(update_bc, field);
             tag.resize();
             AMR_criteria(field, tag);
diff --git a/demos/Weno/heat_weno5_amr.cpp b/demos/Weno/heat_weno5_amr.cpp
index 178c84fc2..1fb32b50b 100644
--- a/demos/Weno/heat_weno5_amr.cpp
+++ b/demos/Weno/heat_weno5_amr.cpp
@@ -1,6 +1,7 @@
 // #include <spdlog/spdlog.h>
 #include <fmt/format.h>
 #include <fmt/ostream.h>
+#include <samurai/print.hpp>
 #include <xtensor/xfixed.hpp>
 #include <xtensor/xmasked_view.hpp>
 
@@ -749,12 +750,12 @@ int main()
 
     for (std::size_t ite = 0; ite < max_ite; ++ite)
     {
-        std::cout << "iteration: " << ite << std::endl;
+        samurai::io::print(samurai::io::root, "iteration: {}\n", ite);
 
         std::size_t ite_adapt = 0;
         while (1)
         {
-            std::cout << "\tmesh adaptation: " << ite_adapt++ << std::endl;
+            samurai::io::print(samurai::io::root, "\tmesh adaptation: {}\n", ite_adapt++);
             samurai::update_ghost(update_bc, field);
             auto tag = samurai::make_scalar_field<int>("tag", mesh);
             AMR_criteria(field, tag);
diff --git a/demos/Weno/weno5_amr.cpp b/demos/Weno/weno5_amr.cpp
index cb7c8bd59..ac31f78bf 100644
--- a/demos/Weno/weno5_amr.cpp
+++ b/demos/Weno/weno5_amr.cpp
@@ -1,6 +1,7 @@
 // #include <spdlog/spdlog.h>
 #include <fmt/format.h>
 #include <fmt/ostream.h>
+#include <samurai/print.hpp>
 #include <xtensor/xfixed.hpp>
 #include <xtensor/xmasked_view.hpp>
 
@@ -464,12 +465,12 @@ int main()
 
     for (std::size_t ite = 0; ite < max_ite; ++ite)
     {
-        std::cout << "iteration: " << ite << std::endl;
+        samurai::io::print("iteration: {}\n", ite);
 
         std::size_t ite_adapt = 0;
         while (1)
         {
-            std::cout << "\tmesh adaptation: " << ite_adapt++ << std::endl;
+            samurai::io::print("\tmesh adaptation: {}\n", ite_adapt++);
             samurai::update_ghost(update_bc, field);
             auto tag = samurai::make_scalar_field<int>("tag", mesh);
             AMR_criteria(field, tag);
diff --git a/demos/highorder/main.cpp b/demos/highorder/main.cpp
index 9c68c95a0..b42ae5e7e 100644
--- a/demos/highorder/main.cpp
+++ b/demos/highorder/main.cpp
@@ -4,6 +4,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -198,7 +199,7 @@ int main(int argc, char* argv[])
             // mesh = {cl, mesh.min_level() + 1, mesh.max_level() + 1};
             mesh = {cl, min_level + i_ref + 1, max_level + i_ref + 1};
         }
-        // std::cout << mesh << std::endl;
+        // samurai::io::print("{}\n", fmt::streamed(mesh));
         // samurai::save("refine_mesh", mesh);
 
         // Equation: -Lap u = f   in [0, 1]^2
@@ -273,29 +274,28 @@ int main(int argc, char* argv[])
         double h = mesh.cell_length(mesh.min_level());
 
         double error = L2_error(u, exact_func);
-        std::cout.precision(2);
-        std::cout << "refinement: " << ite << std::endl;
-        std::cout << "L2-error         : " << std::scientific << error;
+
+        samurai::io::print(samurai::io::root, "refinement: {}\n", ite);
+        samurai::io::print(samurai::io::root, "L2-error         : {:e}", error);
         if (h_coarse != -1)
         {
             double convergence_order = samurai::convergence_order(h, error, h_coarse, error_coarse);
-            std::cout << " (order = " << std::defaultfloat << convergence_order << ")";
+            samurai::io::print(samurai::io::root, " (order = {})", convergence_order);
         }
-        std::cout << std::endl;
+        samurai::io::print(samurai::io::root, "\n");
 
         auto u_recons = samurai::reconstruction(u);
 
         double error_recons = L2_error(u_recons, exact_func);
 
-        std::cout.precision(2);
-        std::cout << "L2-error (recons): " << std::scientific << error_recons;
+        samurai::io::print(samurai::io::root, "L2-error (recons): {:e}", error_recons);
 
         if (h_coarse != -1)
         {
             double convergence_order = samurai::convergence_order(h, error_recons, h_coarse, error_recons_coarse);
-            std::cout << " (order = " << std::defaultfloat << convergence_order << ")";
+            samurai::io::print(samurai::io::root, " (order = {})", convergence_order);
         }
-        std::cout << std::endl;
+        samurai::io::print(samurai::io::root, "\n");
 
         auto error_field = samurai::make_scalar_field<double>("error", mesh);
         samurai::for_each_cell(mesh,
diff --git a/demos/multigrid/main.cpp b/demos/multigrid/main.cpp
index c7eb04e8d..9aa4f7f8d 100644
--- a/demos/multigrid/main.cpp
+++ b/demos/multigrid/main.cpp
@@ -1,12 +1,15 @@
 // Copyright 2018-2025 the samurai's authors
 // SPDX-License-Identifier:  BSD-3-Clause
 
+#include <fmt/core.h>
+#include <fmt/ostream.h>
 #include <iostream>
 #include <samurai/amr/mesh.hpp>
 #include <samurai/box.hpp>
 #include <samurai/field.hpp>
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/schemes/fv.hpp>
 
@@ -166,9 +169,9 @@ int main(int argc, char* argv[])
     }
     else
     {
-        std::cerr << "unknown value for argument --tc" << std::endl;
+        samurai::io::eprint("unknown value for argument --tc\n");
     }
-    std::cout << "Test case: " << test_case_code << std::endl;
+    samurai::io::print("Test case: {}\n", test_case_code);
 
     PetscOptionsGetBool(NULL, NULL, "--save_sol", &save_solution, NULL);
     PetscOptionsGetBool(NULL, NULL, "--save_mesh", &save_mesh, NULL);
@@ -202,11 +205,11 @@ int main(int argc, char* argv[])
 
     if (save_mesh)
     {
-        std::cout << "Saving mesh..." << std::endl;
+        samurai::io::print("Saving mesh...\n");
         samurai::save(path, "mesh", mesh);
     }
 
-    std::cout << "Unknowns: " << (mesh.nb_cells() * n_comp) << std::endl;
+    samurai::io::print("Unknowns: {}\n", (mesh.nb_cells() * n_comp));
 
     //----------------//
     // Create problem //
@@ -243,19 +246,19 @@ int main(int argc, char* argv[])
 
     total_timer.Start();
 
-    std::cout << "Setup solver..." << std::endl;
+    samurai::io::print("Setup solver...\n");
     setup_timer.Start();
     solver.setup();
     setup_timer.Stop();
 
-    std::cout << "Solving..." << std::endl;
+    samurai::io::print("Solving...\n");
     solve_timer.Start();
     solver.solve(source);
     solve_timer.Stop();
 
     total_timer.Stop();
 
-    std::cout << solver.iterations() << " iterations" << std::endl << std::endl;
+    samurai::io::print("{} iterations\n\n", solver.iterations());
 
     solver.destroy_petsc_objects();
 
@@ -263,16 +266,16 @@ int main(int argc, char* argv[])
     //  Print exec times  //
     //--------------------//
 
-    std::cout << "---- Setup ----" << std::endl;
-    std::cout << "CPU time    : " << setup_timer.CPU() << std::endl;
-    std::cout << "Elapsed time: " << setup_timer.Elapsed() << std::endl;
-    std::cout << "---- Solve ----" << std::endl;
-    std::cout << "CPU time    : " << solve_timer.CPU() << std::endl;
-    std::cout << "Elapsed time: " << solve_timer.Elapsed() << std::endl;
-    std::cout << "---- Total ----" << std::endl;
-    std::cout << "CPU time    : " << total_timer.CPU() << std::endl;
-    std::cout << "Elapsed time: " << total_timer.Elapsed() << std::endl;
-    std::cout << std::endl;
+    samurai::io::print("---- Setup ----\n");
+    samurai::io::print("CPU time    : {}\n", fmt::streamed(setup_timer.CPU()));
+    samurai::io::print("Elapsed time: {}\n", fmt::streamed(setup_timer.Elapsed()));
+    samurai::io::print("---- Solve ----\n");
+    samurai::io::print("CPU time    : {}\n", fmt::streamed(solve_timer.CPU()));
+    samurai::io::print("Elapsed time: {}\n", fmt::streamed(solve_timer.Elapsed()));
+    samurai::io::print("---- Total ----\n");
+    samurai::io::print("CPU time    : {}\n", fmt::streamed(total_timer.CPU()));
+    samurai::io::print("Elapsed time: {}\n", fmt::streamed(total_timer.Elapsed()));
+    samurai::io::print("\n");
 
     /*auto right_fluxes = samurai::make_vector_field<double, n_comp, is_soa>("fluxes", mesh);
     samurai::DirectionVector<dim> right = {1, 0};
@@ -292,21 +295,20 @@ int main(int argc, char* argv[])
     if (test_case->solution_is_known())
     {
         double error = L2_error(solution, test_case->solution());
-        std::cout.precision(2);
-        std::cout << "L2-error: " << std::scientific << error << std::endl;
+        samurai::io::print("L2-error: {:.2e}\n", error);
 
         if (test_case_code == "poly")
         {
             // double hidden_constant             = samurai::compute_error_bound_hidden_constant<order>(h, error);
-            // std::cout << "hidden_constant: " << hidden_constant << std::endl;
+            // samurai::io::print("hidden_constant: {}\n", hidden_constant);
             static constexpr std::size_t order = 2;
             double h                           = mesh.cell_length(mesh.min_level());
             double hidden_constant             = 5e-2;
             double theoretical_bound           = samurai::theoretical_error_bound<order>(n_comp * hidden_constant, h);
-            // std::cout << "theoretical_bound: " << theoretical_bound << std::endl;
+            // samurai::io::print("theoretical_bound: {}\n", theoretical_bound);
             if (error > theoretical_bound)
             {
-                std::cerr << "Convergence order failure: the error must be < " << theoretical_bound << "." << std::endl;
+                samurai::io::eprint("Convergence order failure: the error must be < {}.\n", theoretical_bound);
             }
         }
     }
@@ -314,7 +316,7 @@ int main(int argc, char* argv[])
     // Save solution
     if (save_solution)
     {
-        std::cout << "Saving solution..." << std::endl;
+        samurai::io::print("Saving solution...\n");
         samurai::save(path, filename, mesh, solution);
     }
 
diff --git a/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp b/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp
index 1e2603914..8cf6c72d4 100644
--- a/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp
+++ b/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp
@@ -1,5 +1,6 @@
 #pragma once
 #include "SamuraiDM.hpp"
+#include <samurai/print.hpp>
 
 namespace samurai_new
 {
@@ -81,47 +82,47 @@ namespace samurai_new
                     }
                     else
                     {
-                        std::cout << "ERROR: unknown value for argument --smooth" << std::endl << std::endl;
+                        samurai::io::eprint("ERROR: unknown value for argument --smooth\n\n");
                     }
                 }
 
-                std::cout << "Samurai multigrid: " << std::endl;
+                samurai::io::print(samurai::io::root, "Samurai multigrid: \n");
 
-                std::cout << "    smoothers         : ";
+                samurai::io::print(samurai::io::root, "    smoothers         : ");
                 if (smoother == GaussSeidel)
                 {
-                    std::cout << "Gauss-Seidel (pre: lexico., post: antilexico.)";
+                    samurai::io::print(samurai::io::root, "Gauss-Seidel (pre: lexico., post: antilexico.)");
                 }
                 else if (smoother == SymGaussSeidel)
                 {
-                    std::cout << "symmetric Gauss-Seidel";
+                    samurai::io::print(samurai::io::root, "symmetric Gauss-Seidel");
                 }
                 else if (smoother == Petsc)
                 {
-                    std::cout << "petsc options";
+                    samurai::io::print(samurai::io::root, "petsc options");
                 }
-                std::cout << std::endl;
+                samurai::io::print(samurai::io::root, "\n");
 
-                std::cout << "    transfer operators: ";
+                samurai::io::print(samurai::io::root, "    transfer operators: ");
                 if (transfer_ops == TransferOperators::Assembled)
                 {
-                    std::cout << "P assembled, R assembled";
+                    samurai::io::print(samurai::io::root, "P assembled, R assembled");
                 }
                 else if (transfer_ops == TransferOperators::Assembled_PTranspose)
                 {
-                    std::cout << "P assembled, R = P^T";
+                    samurai::io::print(samurai::io::root, "P assembled, R = P^T");
                 }
                 else if (transfer_ops == TransferOperators::MatrixFree_Arrays)
                 {
-                    std::cout << "P mat-free, R mat-free (via double*)";
+                    samurai::io::print(samurai::io::root, "P mat-free, R mat-free (via double*)");
                 }
                 else if (transfer_ops == TransferOperators::MatrixFree_Fields)
                 {
-                    std::cout << "P mat-free, R mat-free (via Fields)";
+                    samurai::io::print(samurai::io::root, "P mat-free, R mat-free (via Fields)");
                 }
-                std::cout << std::endl;
+                samurai::io::print(samurai::io::root, "\n");
 
-                std::cout << "    prediction order  : " << prediction_order << std::endl;
+                samurai::io::print(samurai::io::root, "    prediction order  : {}\n", prediction_order);
 
                 _samuraiDM = new SamuraiDM<Dsctzr>(PETSC_COMM_SELF, *_discretizer, *_mesh, transfer_ops, prediction_order);
                 KSPSetDM(ksp, _samuraiDM->PetscDM());
@@ -147,7 +148,7 @@ namespace samurai_new
                     levels = std::max(static_cast<int>(_mesh->max_level()) - 3, 2);
                     levels = std::min(levels, 8);
                 }
-                std::cout << "    levels            : " << levels << std::endl;
+                samurai::io::print(samurai::io::root, "    levels            : {}\n", levels);
                 PCMGSetLevels(mg, levels, nullptr);
 
                 // All of the following must be called after PCMGSetLevels()
diff --git a/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp b/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp
index 3c56c1458..2f4cd458d 100644
--- a/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp
+++ b/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp
@@ -72,7 +72,7 @@ namespace samurai_new
 
             static PetscErrorCode Coarsen(DM fine_dm, MPI_Comm /*comm*/, DM* coarse_dm)
             {
-                // std::cout << "Coarsen - begin" << std::endl;
+                // samurai::io::print("Coarsen - begin\n");
                 LevelContext<Dsctzr>* fine_ctx;
                 DMShellGetContext(fine_dm, &fine_ctx);
 
@@ -84,10 +84,9 @@ namespace samurai_new
                 DMShellCreate(PetscObjectComm(reinterpret_cast<PetscObject>(fine_dm)), coarse_dm);
                 DefineShellFunctions(*coarse_dm, *coarse_ctx);
 
-                // DMShellCreate(PetscObjectComm((PetscObject)fine_dm),
-                // coarse_dm); DefineShellFunctions(*coarse_dm, *coarse_ctx);
-                // std::cout << "Coarsen (create level " << coarse_ctx->level <<
-                // ")" << std::endl;
+                // DMShellCreate(PetscObjectComm((PetscObject)fine_dm), coarse_dm);
+                // DefineShellFunctions(*coarse_dm, *coarse_ctx);
+                // samurai::io::print("Coarsen (create level {})\n", coarse_ctx->level);
                 return 0;
             }
 
@@ -97,8 +96,7 @@ namespace samurai_new
                 DMShellGetContext(shell, &ctx);
 
                 ctx->assembly().create_matrix(*A);
-                // std::cout << "CreateMatrix - level " << ctx->level <<
-                // std::endl;
+                // samurai::io::print("CreateMatrix - level {}\n", ctx->level);
 
                 MatSetDM(*A, shell); // Why???
                 return 0;
@@ -116,7 +114,7 @@ namespace samurai_new
                 ctx->assembly().assemble_matrix(jac);
 
                 // MatView(jac, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                // std::cout << std::endl;
+                // samurai::io::print("\n");
                 return 0;
             }
 
@@ -147,9 +145,9 @@ namespace samurai_new
                 MatShellGetContext(P, &coarse_ctx);
                 fine_ctx = coarse_ctx->finer;
 
-                // std::cout << "coarse vector x:" << std::endl;
-                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout
-                // << std::endl;
+                // samurai::io::print("coarse vector x:\n");
+                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                // samurai::io::print("\n");
 
                 if (coarse_ctx->transfer_ops == TransferOperators::MatrixFree_Fields)
                 {
@@ -158,9 +156,9 @@ namespace samurai_new
                     Field fine_field = multigrid::prolong(coarse_field, fine_ctx->mesh(), coarse_ctx->prediction_order);
                     samurai::petsc::copy(fine_field, y);
 
-                    // std::cout << "prolongated vector (marche):" << std::endl;
+                    // samurai::io::print("prolongated vector (marche):\n");
                     // VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                    // std::cout << std::endl;
+                    // samurai::io::print("\n");
                 }
                 else
                 {
@@ -172,9 +170,9 @@ namespace samurai_new
                     VecRestoreArrayRead(x, &xarray);
                     VecRestoreArray(y, &yarray);
 
-                    /*std::cout << "prolongated vector (marche):" << std::endl;
-                    VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout
-                    << std::endl;
+                    /*samurai::io::print("prolongated vector (marche):\n");
+                    VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                    samurai::io::print("\n");
 
                     // With matrix
                     std::size_t nf = fine_ctx->mesh().nb_cells();
@@ -191,17 +189,18 @@ namespace samurai_new
                     fine_ctx->mesh(), P2); MatAssemblyBegin(P2,
                     MAT_FINAL_ASSEMBLY); MatAssemblyEnd(P2, MAT_FINAL_ASSEMBLY);
 
-                    std::cout << "prolongated vector with matrix (marche pas):"
-                    << std::endl; Vec y2; VecCreateSeq(PETSC_COMM_SELF, nf,
-                    &y2); MatMult(P2, x, y2); VecView(y2,
-                    PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout <<
-                    std::endl;*/
+                    samurai::io::print("prolongated vector with matrix (marche pas):\n");
+                    Vec y2; VecCreateSeq(PETSC_COMM_SELF, nf, &y2);
+                    MatMult(P2, x, y2);
+                    VecView(y2, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                    samurai::io::print("\n");*/
                 }
 
-                // std::cout << "prolongated vector:" << std::endl;
-                // VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout
-                // << std::endl; PetscReal norm; VecNorm(y, NORM_2, &norm);
-                // std::cout << "prolongated vector norm:" << norm << std::endl;
+                // samurai::io::print("prolongated vector:\n");
+                // VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                // samurai::io::print("\n");
+                // PetscReal norm; VecNorm(y, NORM_2, &norm);
+                // samurai::io::print("prolongated vector norm:{}\n", norm);
 
                 assert(samurai::petsc::check_nan_or_inf(y) && "Nan or Inf after prolongation");
                 return 0;
@@ -231,8 +230,8 @@ namespace samurai_new
                 MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
                 MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);
 
-                // MatView(*P, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout
-                // << std::endl;
+                // MatView(*P, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                // samurai::io::print("\n");
 
                 *scaling = nullptr; // Why???
 
@@ -260,9 +259,9 @@ namespace samurai_new
                 MatShellGetContext(R, &fine_ctx);
                 LevelContext<Dsctzr>* coarse_ctx = fine_ctx->coarser;
 
-                // std::cout << "restriction - fine vector:" << std::endl;
-                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout
-                // << std::endl;
+                // samurai::io::print("restriction - fine vector:\n");
+                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                // samurai::io::print("\n");
 
                 if (coarse_ctx->transfer_ops == TransferOperators::MatrixFree_Fields)
                 {
@@ -282,10 +281,11 @@ namespace samurai_new
                     VecRestoreArray(y, &yarray);
                 }
 
-                // std::cout << "restriction - restricted vector:" << std::endl;
-                // VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout
-                // << std::endl; PetscReal norm; VecNorm(y, NORM_2, &norm);
-                // std::cout << "restricted vector norm:" << norm << std::endl;
+                // samurai::io::print("restriction - restricted vector:\n");
+                // VecView(y, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
+                // samurai::io::print("\n");
+                // PetscReal norm; VecNorm(y, NORM_2, &norm);
+                // samurai::io::print("restricted vector norm:{}\n", norm);
 
                 assert(samurai::petsc::check_nan_or_inf(y) && "Nan or Inf after restriction");
                 return 0;
@@ -320,31 +320,30 @@ namespace samurai_new
                 MatAssemblyEnd(*R, MAT_FINAL_ASSEMBLY);
 
                 // MatView(*R, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                // std::cout << std::endl;
+                // samurai::io::print("\n");
 
                 return 0;
             }
 
             static PetscErrorCode CreateGlobalVector(DM shell, Vec* x)
             {
-                // std::cout << "CreateGlobalVector - begin " << std::endl;
+                // samurai::io::print("CreateGlobalVector - begin\n");
                 LevelContext<Dsctzr>* ctx;
                 DMShellGetContext(shell, &ctx);
                 VecCreateSeq(PETSC_COMM_SELF, static_cast<PetscInt>(ctx->mesh().nb_cells()), x);
                 VecSetDM(*x, shell);
-                // std::cout << "CreateGlobalVector - level " << ctx->level <<
-                // std::endl;
+                // samurai::io::print("CreateGlobalVector - level {}\n", ctx->level);
                 return 0;
             }
 
             static PetscErrorCode CreateLocalVector(DM shell, Vec* x)
             {
-                // std::cout << "CreateLocalVector - begin" << std::endl;
+                // samurai::io::print("CreateLocalVector - begin\n");
                 LevelContext<Dsctzr>* ctx;
                 DMShellGetContext(shell, &ctx);
                 VecCreateSeq(PETSC_COMM_SELF, static_cast<PetscInt>(ctx->mesh().nb_cells()), x);
                 VecSetDM(*x, shell);
-                // std::cout << "CreateLocalVector - end" << std::endl;
+                // samurai::io::print("CreateLocalVector - end\n");
                 return 0;
             }
         };
diff --git a/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp b/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp
index b0d6853a2..bafc2dd9d 100644
--- a/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp
+++ b/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp
@@ -1047,8 +1047,7 @@ namespace samurai_new
                         auto& i, const auto& index)
                             {
                                 auto j = index[0];
-                                std::cout << "avant i = " << i << ", j = " << j
-                        << ":" << coarse_field(level, i, j) << std::endl; if
+                                samurai::io::print("avant i = {}, j = {}:{}\n", i, j, coarse_field(level, i, j)); if
                         (s(0) == 1 && s(1) == 1) // top-right
                                     coarse_field(level, i, j) =
                         fine_field(level+1, 2*i  , 2*j  ); else if (s(0) == -1
@@ -1057,9 +1056,8 @@ namespace samurai_new
                         && s(1) == -1) // bottom-left coarse_field(level, i, j)
                         = fine_field(level+1, 2*i+1, 2*j+1); else if (s(0) == 1
                         && s(1) == -1) // bottom-right coarse_field(level, i, j)
-                        = fine_field(level+1, 2*i  , 2*j+1); std::cout << "apres
-                        i = " << i << ", j = " << j << ":" <<
-                        coarse_field(level, i, j) << std::endl;
+                        = fine_field(level+1, 2*i  , 2*j+1); samurai::io::print("apres i = {}, j = {}:{}\n", i, j, coarse_field(level, i,
+                        j));
                             });
                         }*/
                     }
diff --git a/demos/multigrid/samurai_new/utils.cpp b/demos/multigrid/samurai_new/utils.cpp
index 69846fbe8..ea1a029b1 100644
--- a/demos/multigrid/samurai_new/utils.cpp
+++ b/demos/multigrid/samurai_new/utils.cpp
@@ -1,15 +1,15 @@
 #pragma once
 #include <samurai/algorithm.hpp>
+#include <samurai/print.hpp>
 
 template <class Mesh>
 void print_mesh(Mesh& mesh)
 {
-    std::cout << mesh << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(mesh));
     samurai::for_each_cell(mesh,
                            [](const auto& cell)
                            {
-                               std::cout << "level: " << cell.level << ", cell index: " << cell.index << ", center: " << cell.center(0)
-                                         << std::endl;
+                               samurai::io::print("level: {}, cell index: {}, center: {}\n", cell.level, cell.index, cell.center(0));
                            });
 }
 
diff --git a/demos/p4est/simple_2d.cpp b/demos/p4est/simple_2d.cpp
index e7fea6a03..d14dd2f96 100644
--- a/demos/p4est/simple_2d.cpp
+++ b/demos/p4est/simple_2d.cpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier:  BSD-3-Clause
 
 #include <iostream>
+#include <samurai/print.hpp>
 
 #include <filesystem>
 
@@ -229,14 +230,14 @@ int main(int argc, char* argv[])
     cl[2][{3}].add_interval({2, 4});
 
     samurai::CellArray<dim> mesh_1(cl);
-    std::cout << "nb_cells: " << mesh_1.nb_cells() << "\n";
+    samurai::io::print(samurai::io::root, "nb_cells: {}\n", mesh_1.nb_cells());
 
     tic();
     refine_1(mesh_1, max_level);
     auto duration = toc();
-    std::cout << "Version 1: " << duration << "s" << std::endl;
+    samurai::io::print(samurai::io::root, "Version 1: {}s\n", duration);
     toc();
-    std::cout << "nb_cells: " << mesh_1.nb_cells() << "\n";
+    samurai::io::print(samurai::io::root, "nb_cells: {}\n", mesh_1.nb_cells());
 
     // samurai::CellArray<dim> mesh_2(cl);
     using Config    = samurai::MRConfig<dim>;
@@ -246,12 +247,12 @@ int main(int argc, char* argv[])
     tic();
     refine_2(mesh_2, max_level);
     duration = toc();
-    std::cout << "Version 2: " << duration << "s" << std::endl;
-    std::cout << "nb_cells: " << mesh_2.nb_cells(mesh_id_t::cells) << "\n";
+    samurai::io::print(samurai::io::root, "Version 2: {}s\n", duration);
+    samurai::io::print(samurai::io::root, "nb_cells: {}\n", mesh_2.nb_cells(mesh_id_t::cells));
 
-    std::cout << "Memory used " << std::endl;
+    samurai::io::print(samurai::io::root, "Memory used \n");
     auto mem = samurai::memory_usage(mesh_2, /*verbose*/ true);
-    std::cout << "Total: " << mem << std::endl;
+    samurai::io::print(samurai::io::root, "Total: {}\n", mem);
 
     auto level = samurai::make_scalar_field<std::size_t>("level", mesh_2);
     samurai::for_each_cell(mesh_2,
diff --git a/demos/pablo/bubble_2d.cpp b/demos/pablo/bubble_2d.cpp
index 6a9188ee2..3167f2e90 100644
--- a/demos/pablo/bubble_2d.cpp
+++ b/demos/pablo/bubble_2d.cpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier:  BSD-3-Clause
 
 #include <iostream>
+#include <samurai/print.hpp>
 
 #include <filesystem>
 
@@ -309,7 +310,7 @@ int main(int argc, char* argv[])
         bb_xcenter = bb0_xcenter + aa * xt::cos(omega * t);
         bb_ycenter = bb_ycenter + dt * dy;
 
-        std::cout << fmt::format("iteration -> {} t -> {}", nt++, t) << std::endl;
+        samurai::io::print(samurai::io::root, "iteration -> {} t -> {}\n", nt++, t);
 
         for (std::size_t rep = 0; rep < 10; ++rep)
         {
diff --git a/demos/tutorial/2D_mesh.cpp b/demos/tutorial/2D_mesh.cpp
index 8e08ea98b..4d64d834a 100644
--- a/demos/tutorial/2D_mesh.cpp
+++ b/demos/tutorial/2D_mesh.cpp
@@ -1,7 +1,7 @@
 // Copyright 2018-2025 the samurai's authors
 // SPDX-License-Identifier:  BSD-3-Clause
 
-#include <iostream>
+#include <samurai/print.hpp>
 
 #include <filesystem>
 
@@ -54,7 +54,7 @@ int main(int argc, char* argv[])
 
     const samurai::CellArray<dim> ca{cl};
 
-    std::cout << ca << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(ca));
 
     samurai::save(path, filename, ca);
     samurai::finalize();
diff --git a/demos/tutorial/AMR_1D_Burgers/step_0/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_0/main.cpp
index 925614ecc..4707a0620 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_0/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_0/main.cpp
@@ -6,6 +6,7 @@
 #include <samurai/box.hpp>
 #include <samurai/cell_array.hpp>
 #include <samurai/io/hdf5.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 namespace fs = std::filesystem;
@@ -53,7 +54,7 @@ int main(int argc, char* argv[])
 
     mesh[init_level] = {init_level, box, 0, 1};
 
-    std::cout << mesh << "\n";
+    samurai::io::print("{}\n", fmt::streamed(mesh));
 
     samurai::save(path, filename, mesh);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_1/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_1/main.cpp
index 037641a0e..1f274abde 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_1/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_1/main.cpp
@@ -7,6 +7,7 @@
 #include <samurai/cell_array.hpp>
 #include <samurai/field.hpp>
 #include <samurai/io/hdf5.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 
 #include "init_sol.hpp"
@@ -52,7 +53,7 @@ int main(int argc, char* argv[])
     auto phi = init_sol(mesh);
     /////////////////////////////////
 
-    std::cout << mesh << "\n";
+    samurai::io::print("{}\n", fmt::streamed(mesh));
 
     samurai::save(path, filename, mesh, phi);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_2/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_2/main.cpp
index eac549dc7..0cd490f5e 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_2/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_2/main.cpp
@@ -79,7 +79,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        fmt::print("Iteration = {:4d} Time = {:5.4}\n", nt++, t);
+        samurai::io::print("Iteration = {:4d} Time = {:5.4}\n", nt++, t);
 
         update_sol(dt, phi, phi_np1);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp
index 920d331aa..13ef9ea58 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp
@@ -81,7 +81,7 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        fmt::print("Iteration = {:4d} Time = {:5.4}\n", nt++, t);
+        samurai::io::print("Iteration = {:4d} Time = {:5.4}\n", nt++, t);
 
         update_sol(dt, phi, phi_np1);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp
index 48d4573fd..5a8883770 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp
@@ -90,14 +90,14 @@ int main(int argc, char* argv[])
             t = Tf;
         }
 
-        fmt::print("Iteration = {:4d} Time = {:5.4}\n", nt++, t);
+        samurai::io::print("Iteration = {:4d} Time = {:5.4}\n", nt++, t);
 
         std::size_t i_adapt = 0;
         while (i_adapt < (max_level - min_level + 1))
         {
             auto tag = samurai::make_scalar_field<std::size_t>("tag", mesh);
 
-            fmt::print("adaptation iteration : {:4d}\n", i_adapt++);
+            samurai::io::print("adaptation iteration : {:4d}\n", i_adapt++);
             update_ghost(phi);
             AMR_criterion(phi, tag);
             make_graduation(tag);
diff --git a/demos/tutorial/graduation_case_2.cpp b/demos/tutorial/graduation_case_2.cpp
index 8f6821325..08e73849c 100644
--- a/demos/tutorial/graduation_case_2.cpp
+++ b/demos/tutorial/graduation_case_2.cpp
@@ -12,6 +12,7 @@
 #include <samurai/cell_array.hpp>
 #include <samurai/field.hpp>
 #include <samurai/io/hdf5.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/subset/node.hpp>
 
@@ -69,7 +70,7 @@ int main(int argc, char* argv[])
     std::size_t ite = 0;
     while (true)
     {
-        std::cout << "Iteration for remove intersection: " << ite++ << "\n";
+        samurai::io::print("Iteration for remove intersection: {}\n", ite++);
         auto tag = samurai::make_scalar_field<bool>("tag", ca);
         tag.fill(false);
 
@@ -137,7 +138,7 @@ int main(int argc, char* argv[])
     ite = 0;
     while (true)
     {
-        std::cout << "Iteration for graduation: " << ite++ << "\n";
+        samurai::io::print("Iteration for graduation: {}\n", ite++);
         auto tag = samurai::make_scalar_field<bool>("tag", ca);
         tag.fill(false);
 
diff --git a/demos/tutorial/graduation_case_3.cpp b/demos/tutorial/graduation_case_3.cpp
index 73dcd2904..540042b7c 100644
--- a/demos/tutorial/graduation_case_3.cpp
+++ b/demos/tutorial/graduation_case_3.cpp
@@ -11,6 +11,7 @@
 #include <samurai/cell_array.hpp>
 #include <samurai/field.hpp>
 #include <samurai/io/hdf5.hpp>
+#include <samurai/print.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/subset/node.hpp>
 
@@ -59,7 +60,7 @@ int main(int argc, char* argv[])
     std::size_t ite = 0;
     while (true)
     {
-        std::cout << "Iteration for remove intersection: " << ite++ << "\n";
+        samurai::io::print("Iteration for remove intersection: {}\n", ite++);
 
         auto tag = samurai::make_scalar_field<bool>("tag", ca);
         tag.fill(false);
diff --git a/demos/tutorial/interval.cpp b/demos/tutorial/interval.cpp
index 1d93544ce..496b3b7b6 100644
--- a/demos/tutorial/interval.cpp
+++ b/demos/tutorial/interval.cpp
@@ -1,7 +1,7 @@
 // Copyright 2018-2025 the samurai's authors
 // SPDX-License-Identifier:  BSD-3-Clause
 
-#include <iostream>
+#include <samurai/print.hpp>
 
 #include <samurai/cell_array.hpp>
 #include <samurai/cell_list.hpp>
@@ -20,7 +20,7 @@ int main()
 
     const samurai::CellArray<dim> ca{cl};
 
-    std::cout << ca << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(ca));
 
     constexpr std::size_t start_level = 3;
     const samurai::Box<double, dim> box({-1, -1}, {1, 1});
@@ -28,7 +28,7 @@ int main()
 
     ca_box[start_level] = {start_level, box};
 
-    std::cout << ca_box << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(ca_box));
 
     return 0;
 }
diff --git a/demos/tutorial/reconstruction_1d.cpp b/demos/tutorial/reconstruction_1d.cpp
index 3e0946dc1..ac442695a 100644
--- a/demos/tutorial/reconstruction_1d.cpp
+++ b/demos/tutorial/reconstruction_1d.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/reconstruction.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/uniform_mesh.hpp>
@@ -120,7 +121,7 @@ int main(int argc, char* argv[])
     auto t1            = std::chrono::high_resolution_clock::now();
     auto u_reconstruct = reconstruction(u);
     auto t2            = std::chrono::high_resolution_clock::now();
-    std::cout << "execution time " << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count() << std::endl;
+    samurai::io::print("execution time {}\n", std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count());
 
     auto error = samurai::make_scalar_field<double>("error", u_reconstruct.mesh());
     samurai::for_each_interval(u_reconstruct.mesh(),
diff --git a/demos/tutorial/reconstruction_2d.cpp b/demos/tutorial/reconstruction_2d.cpp
index 8dbd50d88..9e75837ed 100644
--- a/demos/tutorial/reconstruction_2d.cpp
+++ b/demos/tutorial/reconstruction_2d.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/reconstruction.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/uniform_mesh.hpp>
@@ -148,7 +149,7 @@ int main(int argc, char* argv[])
     auto t1            = std::chrono::high_resolution_clock::now();
     auto u_reconstruct = reconstruction(u);
     auto t2            = std::chrono::high_resolution_clock::now();
-    std::cout << "execution time " << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count() << std::endl;
+    samurai::io::print("execution time {}\n", std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count());
 
     auto error = samurai::make_scalar_field<double>("error", u_reconstruct.mesh());
     samurai::for_each_interval(u_reconstruct.mesh(),
diff --git a/demos/tutorial/reconstruction_3d.cpp b/demos/tutorial/reconstruction_3d.cpp
index b2def1b5b..34dc63527 100644
--- a/demos/tutorial/reconstruction_3d.cpp
+++ b/demos/tutorial/reconstruction_3d.cpp
@@ -10,6 +10,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/mr/adapt.hpp>
 #include <samurai/mr/mesh.hpp>
+#include <samurai/print.hpp>
 #include <samurai/reconstruction.hpp>
 #include <samurai/samurai.hpp>
 #include <samurai/uniform_mesh.hpp>
@@ -154,7 +155,7 @@ int main(int argc, char* argv[])
     auto t1            = std::chrono::high_resolution_clock::now();
     auto u_reconstruct = reconstruction(u);
     auto t2            = std::chrono::high_resolution_clock::now();
-    std::cout << "execution time " << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count() << std::endl;
+    samurai::io::print("execution time {}\n", std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count());
 
     auto error = samurai::make_scalar_field<double>("error", u_reconstruct.mesh());
     samurai::for_each_interval(u_reconstruct.mesh(),
diff --git a/demos/tutorial/set_operator.cpp b/demos/tutorial/set_operator.cpp
index 4f46700b1..4664a6812 100644
--- a/demos/tutorial/set_operator.cpp
+++ b/demos/tutorial/set_operator.cpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier:  BSD-3-Clause
 
 #include <iostream>
+#include <samurai/print.hpp>
 
 #include <samurai/cell_array.hpp>
 #include <samurai/cell_list.hpp>
@@ -61,33 +62,33 @@ int main()
 
     ca = {cl, true};
 
-    std::cout << ca << "\n";
+    samurai::io::print(samurai::io::root, "{}\n", fmt::streamed(ca));
 
     auto subset = samurai::intersection(ca[0], ca[1]);
 
     subset(
         [&](const auto& i, auto)
         {
-            std::cout << "intersection found in " << i << std::endl;
+            samurai::io::print(samurai::io::root, "intersection found in {}\n", fmt::streamed(i));
         });
 
     subset.on(0)(
         [&](const auto& i, auto)
         {
-            std::cout << "intersection found in " << i << std::endl;
+            samurai::io::print(samurai::io::root, "intersection found in {}\n", fmt::streamed(i));
         });
 
     subset.on(3)(
         [&](const auto& i, auto)
         {
-            std::cout << "intersection found in " << i << std::endl;
+            samurai::io::print(samurai::io::root, "intersection found in {}\n", fmt::streamed(i));
         });
 
     auto subset_d = samurai::difference(ca[0], ca[1]);
     subset_d(
         [&](const auto& i, auto)
         {
-            std::cout << "difference found in " << i << std::endl;
+            samurai::io::print(samurai::io::root, "difference found in {}\n", fmt::streamed(i));
         });
 
     auto u = samurai::make_scalar_field<double>("u", ca);
@@ -105,11 +106,11 @@ int main()
             u(0, i) = 0.5 * (u(1, 2 * i) + u(1, 2 * i + 1));
         });
 
-    std::cout << u << "\n";
+    samurai::io::print(samurai::io::root, "{}\n", fmt::streamed(u));
 
     subset1.on(0).apply_op(projection(u));
 
-    std::cout << u << "\n";
+    samurai::io::print(samurai::io::root, "{}\n", fmt::streamed(u));
 
     samurai::finalize();
     return 0;
diff --git a/docs/source/reference/snippet/subset_1d.cpp b/docs/source/reference/snippet/subset_1d.cpp
index 8e0bf48ec..3263293fc 100644
--- a/docs/source/reference/snippet/subset_1d.cpp
+++ b/docs/source/reference/snippet/subset_1d.cpp
@@ -1,4 +1,4 @@
-#include <iostream>
+#include <samurai/print.hpp>
 
 #include <samurai/cell_array.hpp>
 #include <samurai/cell_list.hpp>
@@ -24,8 +24,8 @@ int main()
                                u[cell] = cell.indices[0];
                            });
 
-    std::cout << "before projection" << std::endl;
-    std::cout << u << std::endl;
+    samurai::io::print("before projection\n");
+    samurai::io::print("{}\n", fmt::streamed(u));
 
     // Make projection on the intersection
     auto subset = samurai::intersection(ca[0], ca[1]).on(0);
@@ -35,8 +35,8 @@ int main()
             u(0, i) = 0.5 * (u(1, 2 * i) + u(1, 2 * i + 1));
         });
 
-    std::cout << "after projection" << std::endl;
-    std::cout << u << std::endl;
+    samurai::io::print("after projection\n");
+    samurai::io::print("{}\n", fmt::streamed(u));
 
     return 0;
 }
diff --git a/docs/source/tutorial/operator_on_subset.rst b/docs/source/tutorial/operator_on_subset.rst
index 80212ad3f..818838e23 100644
--- a/docs/source/tutorial/operator_on_subset.rst
+++ b/docs/source/tutorial/operator_on_subset.rst
@@ -37,7 +37,7 @@ The output of
 
 .. code-block:: c++
 
-    std::cout << ca << std::endl;
+    samurai::io::print("{}\\n", fmt::streamed(ca));
 
 is
 
@@ -84,7 +84,7 @@ Once again, we use the `operator()` of the subset
 
     set([&](const auto& i, auto)
     {
-        std::cout << "Intersection found in " << i << std::endl;
+        samurai::io::print("Intersection found in {}\\n", i);
     });
 
 This operator takes a lambda function with two parameters: the first one is the interval in the x-direction, whereas the second argument is an array of size `dim - 1` with the coordinates of the other dimensions.
@@ -107,7 +107,7 @@ Then if we want the result on level `0`, we employ
 
     set.on(0)([&](const auto& i, auto)
     {
-        std::cout << "intersection found in " << i << std::endl;
+        samurai::io::print("intersection found in {}\\n", i);
     });
 
 The output is
@@ -123,7 +123,7 @@ It is also possible to compute this intersection on a level which does not actua
 
     set.on(3)([&](const auto& i, auto)
     {
-        std::cout << "intersection found in " << i << std::endl;
+        samurai::io::print("intersection found in {}\\n", i);
     });
 
 The corresponding output is
diff --git a/docs/source/tutorial/snippet/2d_mesh_box.cpp b/docs/source/tutorial/snippet/2d_mesh_box.cpp
index e400f174d..520730f36 100644
--- a/docs/source/tutorial/snippet/2d_mesh_box.cpp
+++ b/docs/source/tutorial/snippet/2d_mesh_box.cpp
@@ -1,4 +1,4 @@
-#include <iostream>
+#include <samurai/print.hpp>
 
 #include <samurai/box.hpp>
 #include <samurai/cell_array.hpp>
@@ -14,7 +14,7 @@ int main()
 
     ca_box[start_level] = {start_level, box};
 
-    std::cout << ca_box << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(ca_box));
 
     samurai::save("2d_mesh_box", ca_box);
     return 0;
diff --git a/docs/source/tutorial/snippet/2d_mesh_representation.cpp b/docs/source/tutorial/snippet/2d_mesh_representation.cpp
index 3a98eac4c..e9e1ca236 100644
--- a/docs/source/tutorial/snippet/2d_mesh_representation.cpp
+++ b/docs/source/tutorial/snippet/2d_mesh_representation.cpp
@@ -1,4 +1,4 @@
-#include <iostream>
+#include <samurai/print.hpp>
 
 #include <samurai/cell_array.hpp>
 #include <samurai/cell_list.hpp>
@@ -31,7 +31,7 @@ int main()
 
     samurai::CellArray<dim> ca{cl};
 
-    std::cout << ca << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(ca));
 
     samurai::save("2d_mesh_representation", ca);
     return 0;
diff --git a/docs/source/tutorial/snippet/interval.cpp b/docs/source/tutorial/snippet/interval.cpp
index 59cab7576..e5d03431a 100644
--- a/docs/source/tutorial/snippet/interval.cpp
+++ b/docs/source/tutorial/snippet/interval.cpp
@@ -1,4 +1,4 @@
-#include <iostream>
+#include <samurai/print.hpp>
 
 #include <samurai/cell_array.hpp>
 #include <samurai/cell_list.hpp>
@@ -16,7 +16,7 @@ int main()
 
     samurai::CellArray<dim> ca{cl};
 
-    std::cout << ca << std::endl;
+    samurai::io::print("{}\n", fmt::streamed(ca));
 
     return 0;
 }
diff --git a/include/samurai/algorithm/graduation.hpp b/include/samurai/algorithm/graduation.hpp
index e15574ddf..1b4de3338 100644
--- a/include/samurai/algorithm/graduation.hpp
+++ b/include/samurai/algorithm/graduation.hpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier:  BSD-3-Clause
 
 #pragma once
+#include <fmt/format.h>
 
 #include <xtensor/xmasked_view.hpp>
 #include <xtensor/xtensor.hpp>
@@ -9,6 +10,7 @@
 #include "../array_of_interval_and_point.hpp"
 #include "../cell_flag.hpp"
 #include "../mesh.hpp"
+#include "../print.hpp"
 #include "../stencil.hpp"
 #include "../subset/node.hpp"
 #include "../subset/utils.hpp"
@@ -301,7 +303,7 @@ namespace samurai
                                          isIntersectionEmpty);
                             break;
                         default:
-                            std::cerr << "Warning: Unsupported number of periodic directions (" << directions.size() << ")." << std::endl;
+                            samurai::io::eprint("Warning: Unsupported number of periodic directions ({}) .\n", directions.size());
                     }
                 }
             }
diff --git a/include/samurai/algorithm/update.hpp b/include/samurai/algorithm/update.hpp
index 196e8fed5..8e16f3033 100644
--- a/include/samurai/algorithm/update.hpp
+++ b/include/samurai/algorithm/update.hpp
@@ -13,6 +13,7 @@
 #include "../field.hpp"
 #include "../numeric/prediction.hpp"
 #include "../numeric/projection.hpp"
+#include "../print.hpp"
 #include "../subset/node.hpp"
 #include "../timers.hpp"
 #include "graduation.hpp"
@@ -180,14 +181,21 @@ namespace samurai
 #ifdef SAMURAI_CHECK_NAN
                                     if (xt::any(xt::isnan(field(children_level, {ii_child, ii_child + 1}, index_child))))
                                     {
-                                        std::cerr << std::endl;
+                                        samurai::io::eprint("\n");
 #ifdef SAMURAI_WITH_MPI
                                         mpi::communicator world;
-                                        std::cerr << "[" << world.rank() << "] ";
+                                        samurai::io::eprint("[{}] ", world.rank());
 #endif
-                                        std::cerr << "NaN found in field(" << children_level << "," << ii_child << "," << index_child
-                                                  << ") during projection of the B.C. into the cell at (" << proj_level << ", " << ii << ", "
-                                                  << index << ")   (dir = " << direction << ", layer = " << layer << ")" << std::endl;
+                                        samurai::io::eprint(
+                                            "NaN found in field({}, {}, {}) during projection of the B.C. into the cell at ({}, {}, {})   (dir = {}, layer = {})\n",
+                                            children_level,
+                                            ii_child,
+                                            fmt::streamed(index_child),
+                                            proj_level,
+                                            ii,
+                                            fmt::streamed(index),
+                                            fmt::streamed(direction),
+                                            layer);
 #ifndef NDEBUG
                                         save(fs::current_path(), "update_ghosts", {true, true}, mesh, field);
 #endif
@@ -217,9 +225,14 @@ namespace samurai
                         // However, it can happen in normal conditions if the domain has a hole, so we don't raise an error in that case.
                         if (mesh.domain().is_box())
                         {
-                            std::cerr << "No children found for the ghost at level " << proj_level << ", i = " << ii
-                                      << ", index = " << index << " during projection of the B.C. into the cell at level " << proj_level
-                                      << ", i=" << i_cell << ", index=" << index << std::endl;
+                            samurai::io::eprint(
+                                "No children found for the ghost at level {}, i = {}, index = {} during projection of the B.C. into the cell at level {}, i = {}, index = {}\n",
+                                proj_level,
+                                ii,
+                                fmt::streamed(index),
+                                proj_level,
+                                fmt::streamed(i_cell),
+                                fmt::streamed(index));
                             save(fs::current_path(), "update_ghosts", {true, true}, mesh, field);
                             assert(false);
                         }
@@ -322,14 +335,18 @@ namespace samurai
 #ifdef SAMURAI_CHECK_NAN
                     if (xt::any(xt::isnan(field(pred_level - 1, i_cell >> 1, index >> 1))))
                     {
-                        std::cerr << std::endl;
+                        samurai::io::eprint("\n");
 #ifdef SAMURAI_WITH_MPI
                         mpi::communicator world;
-                        std::cerr << "[" << world.rank() << "] ";
+                        samurai::io::eprint("[{}] ", world.rank());
 #endif
-                        std::cerr << "NaN found in field(" << (pred_level - 1) << "," << (i_cell >> 1) << "," << (index >> 1)
-                                  << ") during prediction of the B.C. into the cell at (" << pred_level << ", " << ii << ", " << index
-                                  << ") " << std::endl;
+                        samurai::io::eprint("NaN found in field({}, {}, {}) during prediction of the B.C. into the cell at ({}, {}, {})\n",
+                                            (pred_level - 1),
+                                            fmt::streamed(i_cell >> 1),
+                                            fmt::streamed(index >> 1),
+                                            pred_level,
+                                            ii,
+                                            fmt::streamed(index));
 #ifndef NDEBUG
                         samurai::save(fs::current_path(), "update_ghosts", {true, true}, field.mesh(), field);
 #endif
@@ -830,13 +847,13 @@ namespace samurai
                         [&](const auto& i, const auto& index)
                         {
                             // delete cell
-                            std::cout << fmt::format("fall intersection between {} {} on level {} in {} {}",
-                                                     world.rank(),
-                                                     neighbour.rank,
-                                                     level,
-                                                     i,
-                                                     index[0])
-                                      << std::endl;
+                            samurai::io::print(samurai::io::root,
+                                               "fall intersection between {} {} on level {} in {} {}\n",
+                                               world.rank(),
+                                               neighbour.rank,
+                                               level,
+                                               fmt::streamed(i),
+                                               fmt::streamed(index[0]));
                         });
                 }
             }
diff --git a/include/samurai/arguments.hpp b/include/samurai/arguments.hpp
index 0b00c62c3..aa14cf498 100644
--- a/include/samurai/arguments.hpp
+++ b/include/samurai/arguments.hpp
@@ -8,24 +8,24 @@ namespace samurai
 {
     namespace args
     {
-        static bool timers = false;
+        inline bool timers = false;
 #ifdef SAMURAI_WITH_MPI
-        static bool dont_redirect_output = false;
+        inline bool print_root_only = false;
 #endif
-        static int finer_level_flux   = 0;
-        static bool refine_boundary   = false;
-        static bool save_debug_fields = false;
+        inline int finer_level_flux   = 0;
+        inline bool refine_boundary   = false;
+        inline bool save_debug_fields = false;
 
         // MRA arguments
-        static double epsilon    = std::numeric_limits<double>::infinity();
-        static double regularity = std::numeric_limits<double>::infinity();
-        static bool rel_detail   = false;
+        inline double epsilon    = std::numeric_limits<double>::infinity();
+        inline double regularity = std::numeric_limits<double>::infinity();
+        inline bool rel_detail   = false;
     }
 
     inline void read_samurai_arguments(CLI::App& app, int& argc, char**& argv)
     {
 #ifdef SAMURAI_WITH_MPI
-        app.add_flag("--dont-redirect-output", args::dont_redirect_output, "Redirect the output for all ranks different of 0")
+        app.add_flag("--print-root-only", args::print_root_only, "Print on root rank only by default for samurai::io::print/eprint")
             ->capture_default_str()
             ->group("IO");
 #endif
diff --git a/include/samurai/assert_log_trace.hpp b/include/samurai/assert_log_trace.hpp
index 62ede9861..0959539bf 100644
--- a/include/samurai/assert_log_trace.hpp
+++ b/include/samurai/assert_log_trace.hpp
@@ -1,28 +1,29 @@
 #pragma once
 
+#include "print.hpp"
 #include <iostream>
 
-#define SAMURAI_ASSERT(condition, msg)                                                                                \
-    do                                                                                                                \
-    {                                                                                                                 \
-        if (!(condition))                                                                                             \
-        {                                                                                                             \
-            std::cerr << "Assertion failed \nin " << __FILE__ << "\n @line " << __LINE__ << ": " << msg << std::endl; \
-            std::terminate();                                                                                         \
-        }                                                                                                             \
+#define SAMURAI_ASSERT(condition, msg)                                                                                 \
+    do                                                                                                                 \
+    {                                                                                                                  \
+        if (!(condition))                                                                                              \
+        {                                                                                                              \
+            samurai::io::eprint("Assertion failed \nin {} \n @line {}: {}\n", __FILE__, __LINE__, fmt::streamed(msg)); \
+            std::terminate();                                                                                          \
+        }                                                                                                              \
     } while (false)
 
 // #ifdef NDEBUG
-#define SAMURAI_LOG(msg)                                \
-    do                                                  \
-    {                                                   \
-        std::cerr << "SMR::Log:: " << msg << std::endl; \
+#define SAMURAI_LOG(msg)                                            \
+    do                                                              \
+    {                                                               \
+        samurai::io::eprint("SMR::Log:: {}\n", fmt::streamed(msg)); \
     } while (0)
 
-#define SAMURAI_TRACE(msg)                                                        \
-    do                                                                            \
-    {                                                                             \
-        std::cerr << "SMR::Trace[line " << __LINE__ << "] :" << msg << std::endl; \
+#define SAMURAI_TRACE(msg)                                                              \
+    do                                                                                  \
+    {                                                                                   \
+        samurai::io::eprint("SMR::Trace[line {}] :{}\n", __LINE__, fmt::streamed(msg)); \
     } while (0)
 // #else
 // #define MGS_LOG( msg )
diff --git a/include/samurai/bc/apply_field_bc.hpp b/include/samurai/bc/apply_field_bc.hpp
index 37d9689eb..ab48a0378 100644
--- a/include/samurai/bc/apply_field_bc.hpp
+++ b/include/samurai/bc/apply_field_bc.hpp
@@ -2,9 +2,11 @@
 // SPDX-License-Identifier:  BSD-3-Clause
 
 #pragma once
+#include <fmt/format.h>
 
 #include "../boundary.hpp"
 #include "../concepts.hpp"
+#include "../print.hpp"
 #include "polynomial_extrapolation.hpp"
 
 namespace samurai
@@ -44,7 +46,7 @@ namespace samurai
         }
         else
         {
-            std::cerr << "Unknown BC type" << std::endl;
+            samurai::io::eprint("Unknown BC type\n");
             exit(EXIT_FAILURE);
         }
     }
diff --git a/include/samurai/bc/polynomial_extrapolation.hpp b/include/samurai/bc/polynomial_extrapolation.hpp
index 917e33fcb..b64a93681 100644
--- a/include/samurai/bc/polynomial_extrapolation.hpp
+++ b/include/samurai/bc/polynomial_extrapolation.hpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier:  BSD-3-Clause
 
 #pragma once
+#include "../print.hpp"
 #include "bc.hpp"
 
 namespace samurai
@@ -50,8 +51,9 @@ namespace samurai
                     {
                         if (std::isnan(field_value(u, cells[c], field_i)))
                         {
-                            std::cerr << "NaN detected in [" << cells[c]
-                                      << "] when applying polynomial extrapolation to fill the outer ghost [" << ghost << "]." << std::endl;
+                            samurai::io::eprint("NaN detected in [{}] when applying polynomial extrapolation to fill the outer ghost [{}].\n",
+                                                fmt::streamed(cells[c]),
+                                                fmt::streamed(ghost));
                             // save(fs::current_path(), "nan_extrapolation", {true, true}, u.mesh(), u);
                             exit(1);
                         }
diff --git a/include/samurai/box.hpp b/include/samurai/box.hpp
index 8f0dd1c76..b75c4c50e 100644
--- a/include/samurai/box.hpp
+++ b/include/samurai/box.hpp
@@ -3,6 +3,7 @@
 
 #pragma once
 
+#include "print.hpp"
 #include "utils.hpp"
 #include <xtensor/xfixed.hpp>
 #include <xtensor/xio.hpp>
@@ -177,7 +178,7 @@ namespace samurai
         if (d == dim - 1 && box.is_valid())
         {
             boxes.push_back(box);
-            // std::cout << box << std::endl;
+            // samurai::io::print("{}\n", fmt::streamed(box));
         }
 
         difference_impl_rec(box, intersection, d + 1, boxes);
@@ -187,7 +188,7 @@ namespace samurai
         if (d == dim - 1 && box.is_valid() && box != intersection) // The intersection is what we want to remove, so we don't add it
         {
             boxes.push_back(box);
-            // std::cout << box << std::endl;
+            // samurai::io::print("{}\n", fmt::streamed(box));
         }
 
         difference_impl_rec(box, intersection, d + 1, boxes);
@@ -197,7 +198,7 @@ namespace samurai
         if (d == dim - 1 && box.is_valid())
         {
             boxes.push_back(box);
-            // std::cout << box << std::endl;
+            // samurai::io::print("{}\n", fmt::streamed(box));
         }
 
         difference_impl_rec(box, intersection, d + 1, boxes);
@@ -309,9 +310,9 @@ namespace samurai
                 // ... and no tolerance is allowed, we raise an error.
                 if (tol == 0)
                 {
-                    std::cerr << "The box " << box << " cannot be exactly represented with a reasonable cell length. ";
-                    std::cerr << "You can modify the box's dimensions or you can set a tolerance so it can be approximately represented."
-                              << std::endl;
+                    samurai::io::eprint(
+                        "The box {} cannot be exactly represented with a reasonable cell length. You can modify the box's dimensions or you can set a tolerance so it can be approximately represented.\n",
+                        fmt::streamed(box));
                     std::exit(1);
                 }
 
@@ -325,7 +326,7 @@ namespace samurai
                     subdivision_length /= 2;
                     approx_length = xt::ceil(box.length() / subdivision_length) * subdivision_length;
                     error         = xt::abs(approx_length - box.length());
-                    // std::cout << "Approximation error: " << error << std::endl;
+                    // samurai::io::print("Approximation error: {}\n", error);
                 }
             }
             else if (tol > 0)
@@ -340,7 +341,7 @@ namespace samurai
                         subdivision_length *= 2;
                         approx_length = xt::ceil(box.length() / subdivision_length) * subdivision_length;
                         error         = xt::abs(approx_length - box.length());
-                        // std::cout << "Approximation error: " << error << std::endl;
+                        // samurai::io::print("Approximation error: {}\n", error);
                     }
                     subdivision_length /= 2;
                     approx_length = xt::ceil(box.length() / subdivision_length) * subdivision_length;
diff --git a/include/samurai/field.hpp b/include/samurai/field.hpp
index fad5aa5d9..a37a3c899 100644
--- a/include/samurai/field.hpp
+++ b/include/samurai/field.hpp
@@ -12,8 +12,7 @@
 #include <filesystem>
 namespace fs = std::filesystem;
 
-#include <fmt/format.h>
-
+#include <xtensor/xio.hpp>
 #include <xtensor/xtensor.hpp>
 #include <xtensor/xview.hpp>
 
@@ -25,6 +24,7 @@ namespace fs = std::filesystem;
 #include "field_expression.hpp"
 #include "mesh_holder.hpp"
 #include "numeric/gauss_legendre.hpp"
+#include "print.hpp"
 #include "storage/containers.hpp"
 #include "timers.hpp"
 
@@ -137,10 +137,10 @@ namespace samurai
                     {
                         if (std::isnan(this->derived_cast().m_storage.data()[static_cast<std::size_t>(i)]))
                         {
-                            // std::cerr << "READ NaN at level " << level << ", in interval " << interval << std::endl;
+                            // samurai::io::eprint("READ NaN at level {}, in interval {}\n", level, interval);
                             auto ii   = i - interval_tmp.index;
                             auto cell = this->derived_cast().mesh().get_cell(level, static_cast<int>(ii), index...);
-                            std::cerr << "READ NaN in " << cell << std::endl;
+                            samurai::io::eprint("READ NaN in {}\n", fmt::streamed(cell));
                             break;
                         }
                     }
@@ -155,7 +155,7 @@ namespace samurai
             {
                 auto interval_tmp = this->derived_cast().get_interval(level, interval, index);
 
-                // std::cout << "READ OR WRITE: " << level << " " << interval << " " << (... << index) << std::endl;
+                // samurai::io::print("READ OR WRITE: {} {} ...\n", level, interval);
                 return view(m_storage, {interval_tmp.index + interval.start, interval_tmp.index + interval.end, interval.step});
             }
 
@@ -165,7 +165,7 @@ namespace samurai
                                    const xt::xtensor_fixed<interval_value_t, xt::xshape<dim - 1>>& index) const
             {
                 auto interval_tmp = this->derived_cast().get_interval(level, interval, index);
-                // std::cout << "READ: " << level << " " << interval << " " << (... << index) << std::endl;
+                // samurai::io::print("READ: {} {} ...\n", level, interval);
                 auto data = view(m_storage, {interval_tmp.index + interval.start, interval_tmp.index + interval.end, interval.step});
 
 #ifdef SAMURAI_CHECK_NAN
@@ -176,10 +176,10 @@ namespace samurai
                     {
                         if (std::isnan(this->derived_cast().m_storage.data()[static_cast<std::size_t>(i)]))
                         {
-                            // std::cerr << "READ NaN at level " << level << ", in interval " << interval << std::endl;
+                            // samurai::io::eprint("READ NaN at level {}, in interval {}\n", level, interval);
                             auto ii   = i - interval_tmp.index;
                             auto cell = this->derived_cast().mesh().get_cell(level, static_cast<int>(ii), index);
-                            std::cerr << "READ NaN in " << cell << std::endl;
+                            samurai::io::eprint("READ NaN in {}\n", fmt::streamed(cell));
                             break;
                         }
                     }
@@ -254,8 +254,8 @@ namespace samurai
 
                 if (xt::any(xt::isnan(data)))
                 {
-                    // std::cout << data << std::endl;
-                    std::cerr << "READ NaN at level " << level << ", " << interval << std::endl;
+                    // samurai::io::print("{}\n", fmt::streamed(data));
+                    samurai::io::eprint("READ NaN at level {}, {}\n", level, fmt::streamed(interval));
                 }
 #endif
                 return data;
@@ -534,7 +534,8 @@ namespace samurai
             //                   level_field[cell] = cell.level;
             //               });
             // save(fs::current_path(), "mesh_throw", {true, true}, this->mesh(), coords, level_field);
-            (std::cout << ... << index) << std::endl;
+            ((samurai::io::print("{} ", fmt::streamed(index))), ...);
+            samurai::io::print("\n");
             throw std::out_of_range(fmt::format("FIELD ERROR on level {}: try to find interval {}", level, interval));
         }
 
@@ -746,7 +747,7 @@ namespace samurai
 
         if (field1.mesh() != field2.mesh())
         {
-            std::cout << "mesh different" << std::endl;
+            samurai::io::print(samurai::io::root, "mesh different\n");
             return false;
         }
 
@@ -1164,7 +1165,8 @@ namespace samurai
             //                   level_field[cell] = cell.level;
             //               });
             // save(fs::current_path(), "mesh_throw", {true, true}, this->mesh(), coords, level_field);
-            (std::cout << ... << index) << std::endl;
+            ((samurai::io::print("{} ", fmt::streamed(index))), ...);
+            samurai::io::print("\n");
             throw std::out_of_range(fmt::format("FIELD ERROR on level {}: try to find interval {}", level, interval));
         }
 
@@ -1375,7 +1377,7 @@ namespace samurai
 
         if (field1.mesh() != field2.mesh())
         {
-            std::cout << "mesh different" << std::endl;
+            samurai::io::print(samurai::io::root, "mesh different\n");
             return false;
         }
 
diff --git a/include/samurai/io/cgal.hpp b/include/samurai/io/cgal.hpp
index 0b15b03c5..b69b67543 100644
--- a/include/samurai/io/cgal.hpp
+++ b/include/samurai/io/cgal.hpp
@@ -166,7 +166,7 @@ namespace samurai
             CGALMesh cgal_mesh;
             if (!PMP::IO::read_polygon_mesh(input_file, cgal_mesh) || cgal_mesh.is_empty())
             {
-                std::cerr << "Invalid input file." << std::endl;
+                samurai::io::eprint("Invalid input file.\n");
                 return CGALMesh();
             }
 
diff --git a/include/samurai/level_cell_array.hpp b/include/samurai/level_cell_array.hpp
index 48b79c717..2555a8057 100644
--- a/include/samurai/level_cell_array.hpp
+++ b/include/samurai/level_cell_array.hpp
@@ -14,8 +14,8 @@
 #include <boost/serialization/vector.hpp>
 #endif
 
+#include "print.hpp"
 #include <fmt/color.h>
-#include <fmt/format.h>
 
 #include "algorithm.hpp"
 #include "box.hpp"
@@ -639,9 +639,9 @@ namespace samurai
 #ifndef NDEBUG
         if (offset < 0)
         {
-            std::cerr << "Error: Interval not found: level " << m_level << ", i = " << interval << ", index = ";
-            ((std::cerr << index << " "), ...);
-            std::cerr << std::endl;
+            samurai::io::eprint("Error: Interval not found: level {}, i = {}, index = ", m_level, fmt::streamed(interval));
+            ((samurai::io::eprint("{} ", index)), ...);
+            samurai::io::eprint("\n");
         }
 #endif
         return m_cells[0][static_cast<std::size_t>(offset)];
@@ -660,12 +660,12 @@ namespace samurai
 #ifndef NDEBUG
         if (offset < 0)
         {
-            std::cerr << "Error: Interval not found: level " << m_level << ", i = " << interval << ", index =";
+            samurai::io::eprint("Error: Interval not found: level {}, i = {}, index =", m_level, fmt::streamed(interval));
             for (std::size_t d = 0; d < dim - 1; ++d)
             {
-                std::cerr << index[d] << " ";
+                samurai::io::eprint("{} ", index[d]);
             }
-            std::cerr << std::endl;
+            samurai::io::eprint("\n");
         }
 #endif
         return m_cells[0][static_cast<std::size_t>(offset)];
@@ -678,12 +678,12 @@ namespace samurai
 #ifndef NDEBUG
         if (offset < 0)
         {
-            std::cerr << "Error: Interval not found: level " << m_level << ", coord = ";
+            samurai::io::eprint("Error: Interval not found: level {}, coord = ", m_level);
             for (std::size_t d = 0; d < dim; ++d)
             {
-                std::cerr << coord[d] << " ";
+                samurai::io::eprint("{} ", coord[d]);
             }
-            std::cerr << std::endl;
+            samurai::io::eprint("\n");
         }
 #endif
         return m_cells[0][static_cast<std::size_t>(offset)];
@@ -1090,9 +1090,11 @@ namespace samurai
         const double warning_tol = 0.5;
         if (scaling_factor > 0 && xt::any(xt::abs(approx_box.length() - box.length()) >= warning_tol * box.length()))
         {
-            std::cerr << "Warning: the box " << box << " is poorly approximated by " << approx_box << ". ";
-            std::cerr << "This is due to a too large scaling factor (" << scaling_factor
-                      << "). Choose a smaller value for a better approximation." << std::endl;
+            samurai::io::eprint(
+                "Warning: the box {} is poorly approximated by {}. This is due to a too large scaling factor ({}). Choose a smaller value for a better approximation.\n",
+                fmt::streamed(box),
+                fmt::streamed(approx_box),
+                scaling_factor);
         }
         m_scaling_factor = scaling_factor;
 
diff --git a/include/samurai/memory.hpp b/include/samurai/memory.hpp
index 06c0e4e83..a234ff043 100644
--- a/include/samurai/memory.hpp
+++ b/include/samurai/memory.hpp
@@ -5,7 +5,7 @@
 
 #include <numeric>
 
-#include <fmt/format.h>
+#include "print.hpp"
 
 #include "level_cell_array.hpp"
 
@@ -47,7 +47,7 @@ namespace samurai
             std::size_t mem_id = memory_usage(mesh[id]);
             if (verbose)
             {
-                std::cout << fmt::format("Mesh {}: {}", id, mem_id) << std::endl;
+                samurai::io::print("Mesh {}: {}\n", id, mem_id);
             }
             mem += mem_id;
         }
diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp
index 1fa2192f7..57ca272a1 100644
--- a/include/samurai/mesh.hpp
+++ b/include/samurai/mesh.hpp
@@ -6,11 +6,10 @@
 #include <array>
 #include <set>
 
-#include <fmt/format.h>
-
 #include "cell_array.hpp"
 #include "cell_list.hpp"
 #include "domain_builder.hpp"
+#include "print.hpp"
 #include "static_algorithm.hpp"
 #include "stencil.hpp"
 #include "subset/node.hpp"
@@ -292,7 +291,7 @@ namespace samurai
         m_periodic.fill(false);
 
 #ifdef SAMURAI_WITH_MPI
-        std::cerr << "MPI is not implemented with DomainBuilder." << std::endl;
+        samurai::io::eprint("MPI is not implemented with DomainBuilder.\n");
         std::exit(1);
         // partition_mesh(start_level, b);
         //  load_balancing();
@@ -482,10 +481,12 @@ namespace samurai
             {
                 if (box.min_length() < 2 * largest_cell_length * config::max_stencil_width)
                 {
-                    std::cerr << "The hole " << box << " is too small to apply the BC at level " << min_level_bc
-                              << " with the given scaling factor. We need to be able to construct " << (2 * config::max_stencil_width)
-                              << " ghosts in each direction inside the hole." << std::endl;
-                    std::cerr << "Please choose a smaller scaling factor or enlarge the hole." << std::endl;
+                    samurai::io::eprint(
+                        "The hole {} is too small to apply the BC at level {} with the given scaling factor. We need to be able to construct {} ghosts in each direction inside the hole.\n",
+                        fmt::streamed(box),
+                        min_level_bc,
+                        (2 * config::max_stencil_width));
+                    samurai::io::eprint("Please choose a smaller scaling factor or enlarge the hole.\n");
                     std::exit(1);
                 }
             }
@@ -1151,7 +1152,7 @@ namespace samurai
             world.barrier();
             if (rank == 0)
             {
-                std::cout << "---------------- k = " << k << " ----------------" << std::endl;
+                samurai::io::print("---------------- k = {} ----------------\n", k);
             }
             mpi::all_gather(world, load, loads);
 
@@ -1181,7 +1182,7 @@ namespace samurai
 
             load_transfer(load_fluxes);
 
-            std::cout << rank << ": load = " << load << ", load_np1 = " << load_np1 << std::endl;
+            samurai::io::print("{}: load = {}, load_np1 = {}\n", rank, load, load_np1);
 
             load = static_cast<std::size_t>(load_np1);
         }
@@ -1193,7 +1194,7 @@ namespace samurai
     {
 #ifdef SAMURAI_WITH_MPI
         mpi::communicator world;
-        std::cout << world.rank() << ": ";
+        samurai::io::print("{}: ", world.rank());
         for (std::size_t i_rank = 0; i_rank < m_mpi_neighbourhood.size(); ++i_rank)
         {
             auto neighbour = m_mpi_neighbourhood[i_rank];
@@ -1203,9 +1204,9 @@ namespace samurai
             else if (load_fluxes[i_rank] > 0) // must receive load from the neighbour
             {
             }
-            std::cout << "--> " << neighbour.rank << ": " << load_fluxes[i_rank] << ", ";
+            samurai::io::print("--> {}: {}, ", neighbour.rank, load_fluxes[i_rank]);
         }
-        std::cout << std::endl;
+        samurai::io::print("\n");
 #endif
     }
 
diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp
index 958c14a7f..0e234ccc7 100644
--- a/include/samurai/mr/adapt.hpp
+++ b/include/samurai/mr/adapt.hpp
@@ -154,7 +154,7 @@ namespace samurai
         cfg.parse_args();
         for (std::size_t i = 0; i < max_level - min_level; ++i)
         {
-            // std::cout << "MR mesh adaptation " << i << std::endl;
+            // samurai::io::print("MR mesh adaptation {}\n", i);
             m_detail.resize();
             m_detail.fill(0);
             m_tag.resize();
@@ -418,7 +418,7 @@ namespace samurai
         requires(!IsField<Prediction_fn>) && (IsField<TFields> && ...)
     auto make_MRAdapt(Prediction_fn&& prediction_fn, TFields&... fields)
     {
-        std::cout << "Use custom prediction function for MRAdapt" << std::endl;
+        samurai::io::print(samurai::io::root, "Use custom prediction function for MRAdapt\n");
         return Adapt<false, Prediction_fn, TFields...>(std::forward<Prediction_fn>(prediction_fn), fields...);
     }
 
diff --git a/include/samurai/mr/operators.hpp b/include/samurai/mr/operators.hpp
index a533680ec..6a42b0e70 100644
--- a/include/samurai/mr/operators.hpp
+++ b/include/samurai/mr/operators.hpp
@@ -11,6 +11,7 @@
 #include "../field.hpp"
 #include "../numeric/prediction.hpp"
 #include "../operators_base.hpp"
+#include "../print.hpp"
 #include "../utils.hpp"
 
 namespace samurai
@@ -288,7 +289,7 @@ namespace samurai
                     {
                         if (std::isnan(qs_i(ii)) || std::isnan(qs_j(ii)) || std::isnan(qs_ij(ii)))
                         {
-                            std::cerr << "NaN detected during the computation of details." << std::endl;
+                            samurai::io::eprint("NaN detected during the computation of details.\n");
                             save(fs::current_path(), "check_nan", {true, true}, field.mesh(), field);
                             break;
                         }
@@ -298,7 +299,7 @@ namespace samurai
                 {
                     if (xt::any(xt::isnan(qs_ij)))
                     {
-                        std::cerr << "NaN detected during the computation of details." << std::endl;
+                        samurai::io::eprint("NaN detected during the computation of details.\n");
                         save(fs::current_path(), "check_nan", {true, true}, field.mesh(), field);
                     }
                 }
diff --git a/include/samurai/numeric/prediction.hpp b/include/samurai/numeric/prediction.hpp
index d95209fc9..b46da1516 100644
--- a/include/samurai/numeric/prediction.hpp
+++ b/include/samurai/numeric/prediction.hpp
@@ -6,6 +6,7 @@
 #include <array>
 #include <type_traits>
 
+#include "../print.hpp"
 #include <xtensor/xadapt.hpp>
 #include <xtensor/xview.hpp>
 
@@ -576,8 +577,7 @@ namespace samurai
 #ifdef SAMURAI_CHECK_NAN
         if (xt::any(xt::isnan(qs_ij)))
         {
-            std::cerr << "NaN detected in the prediction stencil (Qs_ij)." << std::endl;
-            std::cerr << qs_ij << std::endl;
+            samurai::io::eprint("NaN detected in the prediction stencil (Qs_ij).\n{}\n", fmt::streamed(qs_ij));
             save(fs::current_path(), "check_nan", {true, true}, src.mesh(), src);
         }
 #endif
diff --git a/include/samurai/petsc/block_assembly.hpp b/include/samurai/petsc/block_assembly.hpp
index ec0311c06..9e6f81840 100644
--- a/include/samurai/petsc/block_assembly.hpp
+++ b/include/samurai/petsc/block_assembly.hpp
@@ -1,7 +1,9 @@
 #pragma once
+#include "../print.hpp"
 #include "../schemes/block_operator.hpp"
 #include "utils.hpp"
 #include "zero_block_assembly.hpp"
+#include <fmt/format.h>
 
 namespace samurai
 {
@@ -121,8 +123,11 @@ namespace samurai
                                 }
                                 else if (op.matrix_rows() != block_rows)
                                 {
-                                    std::cerr << "Assembly failure: incompatible number of rows of block (" << row << ", " << col
-                                              << "): " << op.matrix_rows() << " (expected " << block_rows << ")" << std::endl;
+                                    samurai::io::eprint("Assembly failure: incompatible number of rows of block ({}, {}): {} (expected {})\n",
+                                                        row,
+                                                        col,
+                                                        op.matrix_rows(),
+                                                        block_rows);
                                     exit(EXIT_FAILURE);
                                 }
                             }
@@ -154,8 +159,12 @@ namespace samurai
                                 }
                                 else if (op.matrix_cols() != block_cols)
                                 {
-                                    std::cerr << "Assembly failure: incompatible number of columns of block (" << row << ", " << col
-                                              << "): " << op.matrix_cols() << " (expected " << block_cols << ")" << std::endl;
+                                    samurai::io::eprint(
+                                        "Assembly failure: incompatible number of columns of block ({}, {}): {} (expected {})\n",
+                                        row,
+                                        col,
+                                        op.matrix_cols(),
+                                        block_cols);
                                     exit(EXIT_FAILURE);
                                 }
                             }
@@ -192,8 +201,11 @@ namespace samurai
                                              }
                                              else
                                              {
-                                                 std::cerr << "unknown " << i << " is not compatible with the scheme (" << row << ", "
-                                                           << col << ") (named '" << op.name() << "')" << std::endl;
+                                                 samurai::io::eprint("unknown {} is not compatible with the scheme ({}, {}) (named '{}')\n",
+                                                                     i,
+                                                                     row,
+                                                                     col,
+                                                                     op.name());
                                                  assert(false);
                                                  exit(EXIT_FAILURE);
                                              }
@@ -241,8 +253,10 @@ namespace samurai
                     {
                         if (row == col && op.matrix_rows() != op.matrix_cols())
                         {
-                            std::cerr << "Function 'create_vector()' is ambiguous in this context, because the block (" << row << ", " << col
-                                      << ") is not square. Use 'create_applicable_vector()' or 'create_rhs_vector()' instead." << std::endl;
+                            samurai::io::eprint(
+                                "Function 'create_vector()' is ambiguous in this context, because the block ({}, {}) is not square. Use 'create_applicable_vector()' or 'create_rhs_vector()' instead.\n",
+                                row,
+                                col);
                             exit(EXIT_FAILURE);
                         }
                     });
@@ -293,7 +307,7 @@ namespace samurai
                 for_each_assembly_op(
                     [&](auto& op, auto row, auto col)
                     {
-                        // std::cout << "create_matrix (" << row << ", " << col << ")" << std::endl;
+                        // samurai::io::print("create_matrix ({}, {})\n", row, col);
                         op.create_matrix(block(row, col));
                     });
                 MatCreateNest(PETSC_COMM_SELF, rows, PETSC_IGNORE, cols, PETSC_IGNORE, m_blocks.data(), &A);
@@ -304,7 +318,7 @@ namespace samurai
                 for_each_assembly_op(
                     [&](auto& op, auto row, auto col)
                     {
-                        // std::cout << "assemble_matrix (" << row << ", " << col << ") '" << op.name() << "'" << std::endl;
+                        // samurai::io::print("assemble_matrix ({}, {}) '{}'\n", row, col, op.name());
                         op.assemble_matrix(block(row, col));
                     });
                 if (final_assembly)
@@ -357,7 +371,7 @@ namespace samurai
                     {
                         if (op.include_bc())
                         {
-                            // std::cout << "enforce_bc (" << row << ", " << col << ") on b[" << row << "]" << std::endl;
+                            // samurai::io::print("enforce_bc ({}, {}) on b[{}]\n", row, col, row);
                             Vec b_block;
                             VecNestGetSubVec(b, static_cast<PetscInt>(row), &b_block);
                             op.enforce_bc(b_block);
@@ -633,7 +647,7 @@ namespace samurai
                         {
                             op.set_current_insert_mode(insert_mode);
                         }
-                        // std::cout << "assemble_scheme (" << row << ", " << col << ") '" << op.name() << "'" << std::endl;
+                        // samurai::io::print("assemble_scheme ({}, {}) '{}'\n", row, col, op.name());
                         op.assemble_scheme(A);
                         insert_mode = op.current_insert_mode();
                     });
@@ -778,7 +792,7 @@ namespace samurai
                     {
                         if (op.include_bc())
                         {
-                            // std::cout << "enforce_bc (" << row << ", " << col << ") on b[" << row << "]" << std::endl;
+                            // samurai::io::print("enforce_bc ({}, {}) on b[{}]\n", row, col, row);
                             op.enforce_bc(b);
                         }
                     });
diff --git a/include/samurai/petsc/fv/FV_scheme_assembly.hpp b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
index abe9ab2b9..bf2dbd31d 100644
--- a/include/samurai/petsc/fv/FV_scheme_assembly.hpp
+++ b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
@@ -2,6 +2,7 @@
 #include "../../algorithm/update.hpp"
 #include "../../boundary.hpp"
 #include "../../numeric/prediction.hpp"
+#include "../../print.hpp"
 #include "../../schemes/fv/FV_scheme.hpp"
 #include "../../schemes/fv/scheme_operators.hpp"
 #include "../matrix_assembly.hpp"
@@ -90,12 +91,12 @@ namespace samurai
             {
                 if (!m_unknown)
                 {
-                    std::cerr << "Undefined unknown for operator " << scheme().name() << ". Assembly initialization failed!" << std::endl;
+                    samurai::io::eprint("Undefined unknown for operator {}. Assembly initialization failed!\n", this->name());
                     assert(false);
                     exit(EXIT_FAILURE);
                 }
                 m_n_cells = mesh().nb_cells();
-                // std::cout << "reset " << this->name() << ", rows = " << matrix_rows() << std::endl;
+                // samurai::io::print("reset {}, rows = {}\n", this->name(), matrix_rows());
                 m_is_row_empty.resize(static_cast<std::size_t>(matrix_rows()));
                 std::fill(m_is_row_empty.begin(), m_is_row_empty.end(), true);
 
@@ -280,8 +281,10 @@ namespace samurai
 
                     if (cfg_t::stencil_size > 1 && unknown().get_bc().empty())
                     {
-                        std::cerr << "Failure to assemble the boundary conditions in the operator '" << this->name()
-                                  << "': no boundary condition attached to the field '" << unknown().name() << "'." << std::endl;
+                        samurai::io::eprint(
+                            "Failure to assemble the boundary conditions in the operator '{}': no boundary condition attached to the field '{}'.\n",
+                            this->name(),
+                            unknown().name());
                         assert(false);
                         continue;
                     }
@@ -329,8 +332,8 @@ namespace samurai
                                 }
                                 else
                                 {
-                                    std::cerr << "Unknown boundary condition type. Only Dirichlet and Neumann are implemented at the moment."
-                                              << std::endl;
+                                    samurai::io::eprint(
+                                        "Unknown boundary condition type. Only Dirichlet and Neumann are implemented at the moment.\n");
                                 }
                             }
                         }
@@ -427,7 +430,7 @@ namespace samurai
 
             void assemble_boundary_conditions(Mat& A) override
             {
-                // std::cout << "assemble_boundary_conditions of " << this->name() << std::endl;
+                // samurai::io::print("assemble_boundary_conditions of {}\n", this->name());
                 if (current_insert_mode() == ADD_VALUES)
                 {
                     // Must flush to use INSERT_VALUES instead of ADD_VALUES
@@ -509,15 +512,18 @@ namespace samurai
 
             virtual void enforce_bc(Vec& b) const
             {
-                // std::cout << "enforce_bc of " << this->name() << std::endl;
+                // samurai::io::print("enforce_bc of {}\n", this->name());
                 if (!this->is_block())
                 {
                     PetscInt b_rows;
                     VecGetSize(b, &b_rows);
                     if (b_rows != this->matrix_cols())
                     {
-                        std::cerr << "Operator '" << this->name() << "': the number of rows in vector (" << b_rows
-                                  << ") does not equal the number of columns of the matrix (" << this->matrix_cols() << ")" << std::endl;
+                        samurai::io::eprint(
+                            "Operator '{}': the number of rows in vector ({}) does not equal the number of columns of the matrix ({})\n",
+                            this->name(),
+                            b_rows,
+                            this->matrix_cols());
                         assert(false);
                         return;
                     }
@@ -594,7 +600,7 @@ namespace samurai
                     MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
                     set_current_insert_mode(INSERT_VALUES);
                 }
-                // std::cout << "insert_value_on_diag_for_useless_ghosts of " << this->name() << std::endl;
+                // samurai::io::print("insert_value_on_diag_for_useless_ghosts of {}\n", this->name());
                 for (std::size_t i = 0; i < m_is_row_empty.size(); i++)
                 {
                     if (m_is_row_empty[i])
@@ -606,9 +612,12 @@ namespace samurai
                                                  INSERT_VALUES);
                         if (error)
                         {
-                            std::cerr << scheme().name() << ": failure to insert diagonal coefficient at ("
-                                      << m_row_shift + static_cast<PetscInt>(i) << ", " << m_col_shift + static_cast<PetscInt>(i)
-                                      << "), i.e. (" << i << ", " << i << ") in the block." << std::endl;
+                            samurai::io::eprint("{}: failure to insert diagonal coefficient at ({}, {}), i.e. ({}, {}) in the block.\n",
+                                                scheme().name(),
+                                                m_row_shift + static_cast<PetscInt>(i),
+                                                m_col_shift + static_cast<PetscInt>(i),
+                                                i,
+                                                i);
                             assert(false);
                             exit(EXIT_FAILURE);
                         }
@@ -753,8 +762,10 @@ namespace samurai
                                                              current_insert_mode());
                                     if (error)
                                     {
-                                        std::cerr << scheme().name() << ": failure to insert projection coefficient at (" << ghost_index
-                                                  << ", " << m_col_shift + col_index(children[i], field_i) << ")." << std::endl;
+                                        samurai::io::eprint("{}: failure to insert projection coefficient at ({}, {}).\n",
+                                                            scheme().name(),
+                                                            ghost_index,
+                                                            m_col_shift + col_index(children[i], field_i));
                                         assert(false);
                                         exit(EXIT_FAILURE);
                                     }
diff --git a/include/samurai/petsc/fv/cell_based_scheme_assembly.hpp b/include/samurai/petsc/fv/cell_based_scheme_assembly.hpp
index 929d313e2..c988ac5ea 100644
--- a/include/samurai/petsc/fv/cell_based_scheme_assembly.hpp
+++ b/include/samurai/petsc/fv/cell_based_scheme_assembly.hpp
@@ -159,7 +159,7 @@ namespace samurai
 
             void assemble_scheme(Mat& A) override
             {
-                // std::cout << "assemble_scheme() of " << this->name() << std::endl;
+                // samurai::io::print("assemble_scheme() of {}\n", this->name());
 
                 if (this->current_insert_mode() == INSERT_VALUES)
                 {
@@ -174,9 +174,9 @@ namespace samurai
                     unknown(),
                     [&](const auto& cells, const auto& coeffs)
                     {
-                        // std::cout << "coeffs: " << std::endl;
+                        // samurai::io::print("coeffs:\n");
                         // for (std::size_t i=0; i<stencil_size; i++)
-                        //     std::cout << i << ": " << coeffs[i] << std::endl;
+                        //     samurai::io::print("{}: {}\n", i, coeffs[i]);
 
                         // Global rows and columns
                         std::array<PetscInt, cfg_t::stencil_size * output_n_comp> rows;
diff --git a/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp b/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp
index fb6434efa..b817737b0 100644
--- a/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp
+++ b/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp
@@ -134,7 +134,7 @@ namespace samurai
 
             void assemble_scheme(Mat& A) override
             {
-                // std::cout << "assemble_scheme() of " << this->name() << std::endl;
+                // samurai::io::print("assemble_scheme() of {}\n", this->name());
 
                 if (this->current_insert_mode() == INSERT_VALUES)
                 {
diff --git a/include/samurai/petsc/fv/operator_sum_assembly.hpp b/include/samurai/petsc/fv/operator_sum_assembly.hpp
index 78fbc4272..c5d15b2b4 100644
--- a/include/samurai/petsc/fv/operator_sum_assembly.hpp
+++ b/include/samurai/petsc/fv/operator_sum_assembly.hpp
@@ -1,6 +1,8 @@
 #pragma once
+#include "../../print.hpp"
 #include "../../schemes/fv/scheme_operators.hpp"
 #include "../matrix_assembly.hpp"
+#include <fmt/format.h>
 
 namespace samurai
 {
@@ -130,10 +132,12 @@ namespace samurai
                          {
                              if (op.matrix_rows() != rows)
                              {
-                                 std::cerr << "Invalid '+' operation: all schemes must generate the same number of matrix rows." << std::endl;
-                                 std::cerr << "                       '" << std::get<0>(m_assembly_ops).name()
-                                           << "': " << std::get<0>(m_assembly_ops).matrix_rows() << ", " << op.name() << ": "
-                                           << op.matrix_rows() << std::endl;
+                                 samurai::io::eprint("Invalid '+' operation: all schemes must generate the same number of matrix rows.\n");
+                                 samurai::io::eprint("                       '{}' : {}, {}: {}\n",
+                                                     std::get<0>(m_assembly_ops).name(),
+                                                     std::get<0>(m_assembly_ops).matrix_rows(),
+                                                     op.name(),
+                                                     op.matrix_rows());
                                  assert(false);
                              }
                          });
@@ -148,11 +152,12 @@ namespace samurai
                          {
                              if (op.matrix_cols() != cols)
                              {
-                                 std::cerr << "Invalid '+' operation: all schemes must generate the same number of matrix columns."
-                                           << std::endl;
-                                 std::cerr << "                       '" << std::get<0>(m_assembly_ops).name()
-                                           << "': " << std::get<0>(m_assembly_ops).matrix_cols() << ", " << op.name() << ": "
-                                           << op.matrix_cols() << std::endl;
+                                 samurai::io::eprint("Invalid '+' operation: all schemes must generate the same number of matrix columns.\n");
+                                 samurai::io::eprint("                       '{}' : {}, {}: {}\n",
+                                                     std::get<0>(m_assembly_ops).name(),
+                                                     std::get<0>(m_assembly_ops).matrix_cols(),
+                                                     op.name(),
+                                                     op.matrix_cols());
                                  assert(false);
                              }
                          });
@@ -188,7 +193,7 @@ namespace samurai
                 for_each(m_assembly_ops,
                          [&](auto& op)
                          {
-                             // std::cout << "Assembly of " << op.name() << std::endl;
+                             // samurai::io::print("Assembly of {}\n", op.name());
                              op.assemble_scheme(A);
                              set_current_insert_mode(op.current_insert_mode());
                          });
diff --git a/include/samurai/petsc/linear_block_solver.hpp b/include/samurai/petsc/linear_block_solver.hpp
index bcbe16f0c..6d97353a2 100644
--- a/include/samurai/petsc/linear_block_solver.hpp
+++ b/include/samurai/petsc/linear_block_solver.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include "../print.hpp"
 #include "block_assembly.hpp"
 #include "linear_solver.hpp"
 
@@ -71,7 +72,7 @@ namespace samurai
                 {
                     if (m_A == nullptr)
                     {
-                        std::cerr << "The matrix must be assemble before calling set_pc_fieldsplit()." << std::endl;
+                        samurai::io::eprint("The matrix must be assemble before calling set_pc_fieldsplit().\n");
                         exit(EXIT_FAILURE);
                     }
                     IS IS_fields[cols];
@@ -115,7 +116,7 @@ namespace samurai
                 }
 
                 this->assemble_matrix();
-                // MatView(m_A, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout << std::endl;
+                // MatView(m_A, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); samurai::io::print("\n");
                 KSPSetOperators(m_ksp, m_A, m_A);
 
                 PC pc;
diff --git a/include/samurai/petsc/linear_solver.hpp b/include/samurai/petsc/linear_solver.hpp
index 28afc89e9..5b889998f 100644
--- a/include/samurai/petsc/linear_solver.hpp
+++ b/include/samurai/petsc/linear_solver.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 // #define ENABLE_MG
+#include "../print.hpp"
 #include "fv/cell_based_scheme_assembly.hpp"
 #include "fv/flux_based_scheme_assembly.hpp"
 #include "fv/operator_sum_assembly.hpp"
@@ -129,8 +130,8 @@ namespace samurai
                 {
                     if (assembly().undefined_unknown())
                     {
-                        std::cerr << "Undefined unknown(s) for this linear system. Please set the unknowns using the instruction '[solver].set_unknown(u);' or '[solver].set_unknowns(u1, u2...);'."
-                                  << std::endl;
+                        samurai::io::eprint(
+                            "Undefined unknown(s) for this linear system. Please set the unknowns using the instruction '[solver].set_unknown(u);' or '[solver].set_unknowns(u1, u2...).\n");
                         assert(false && "Undefined unknown(s)");
                         exit(EXIT_FAILURE);
                     }
@@ -161,7 +162,7 @@ namespace samurai
 
                 if (err != PETSC_SUCCESS)
                 {
-                    std::cerr << "The setup of the solver failed!" << std::endl;
+                    samurai::io::eprint("The setup of the solver failed!\n");
                     assert(false && "Failed solver setup");
                     exit(EXIT_FAILURE);
                 }
@@ -181,7 +182,7 @@ namespace samurai
                 // assembly().enforce_projection_prediction(b);
                 // // Set to zero the right-hand side of the useless ghosts' equations
                 // assembly().set_0_for_useless_ghosts(b);
-                // VecView(b, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout << std::endl;
+                // VecView(b, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); samurai::io::print("\n");
                 // assert(check_nan_or_inf(b));
                 times::timers.stop("system solve");
 
@@ -204,14 +205,14 @@ namespace samurai
                     using namespace std::string_literals;
                     const char* reason_text;
                     KSPGetConvergedReasonString(m_ksp, &reason_text);
-                    std::cerr << "Divergence of the solver ("s + reason_text + ")" << std::endl;
+                    samurai::io::eprint("Divergence of the solver ({})\n", reason_text);
                     // VecView(b, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                    // std::cout << std::endl;
+                    // samurai::io::print("\n");
                     // assert(check_nan_or_inf(b));
                     assert(false && "Divergence of the solver");
                     exit(EXIT_FAILURE);
                 }
-                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout << std::endl;
+                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); samurai::io::print("\n");
             }
 
           public:
@@ -290,9 +291,7 @@ namespace samurai
                 {
                     if constexpr (Mesh::dim > 2)
                     {
-                        std::cerr << "Samurai Multigrid is not implemented for "
-                                     "dim > 2."
-                                  << std::endl;
+                        samurai::io::eprint("Samurai Multigrid is not implemented for dim > 2.\n");
                         assert(false);
                         exit(EXIT_FAILURE);
                     }
@@ -325,8 +324,8 @@ namespace samurai
                 }
                 if (assembly().undefined_unknown())
                 {
-                    std::cerr << "Undefined unknown for this linear system. Please set the unknown using the instruction '[solver].set_unknown(u);'."
-                              << std::endl;
+                    samurai::io::eprint(
+                        "Undefined unknown for this linear system. Please set the unknown using the instruction '[solver].set_unknown(u);'.\n");
                     assert(false && "Undefined unknown");
                     exit(EXIT_FAILURE);
                 }
diff --git a/include/samurai/petsc/matrix_assembly.hpp b/include/samurai/petsc/matrix_assembly.hpp
index 2366da65d..b69587281 100644
--- a/include/samurai/petsc/matrix_assembly.hpp
+++ b/include/samurai/petsc/matrix_assembly.hpp
@@ -194,7 +194,7 @@ namespace samurai
 
                 // for (std::size_t row = 0; row < nnz.size(); ++row)
                 // {
-                //     std::cout << "nnz[" << row << "] = " << nnz[row] << std::endl;
+                //     samurai::io::print("nnz[{}] = {}\n", row, nnz[row]);
                 // }
                 if (!m_is_block)
                 {
@@ -246,7 +246,7 @@ namespace samurai
 
             virtual ~MatrixAssembly()
             {
-                // std::cout << "Destruction of '" << name() << "'" << std::endl;
+                // samurai::io::print("Destruction of '{}'\n", name());
                 m_is_deleted = true;
             }
 
diff --git a/include/samurai/petsc/nonlinear_local_solvers.hpp b/include/samurai/petsc/nonlinear_local_solvers.hpp
index 6bd6201ff..dcfa2cb91 100644
--- a/include/samurai/petsc/nonlinear_local_solvers.hpp
+++ b/include/samurai/petsc/nonlinear_local_solvers.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../print.hpp"
 #include "fv/cell_based_scheme_assembly.hpp"
 #include "fv/flux_based_scheme_assembly.hpp"
 #include "fv/operator_sum_assembly.hpp"
@@ -50,15 +51,14 @@ namespace samurai
             {
                 if (!m_scheme.scheme_definition().local_scheme_function)
                 {
-                    std::cerr << "The scheme function 'local_scheme_function' of operator '" << scheme.name()
-                              << "' has not been implemented." << std::endl;
+                    samurai::io::eprint("The scheme function 'local_scheme_function' of operator '{}' has not been implemented.\n",
+                                        scheme.name());
                     assert(false && "Undefined 'local_scheme_function'");
                     exit(EXIT_FAILURE);
                 }
                 if (!m_scheme.scheme_definition().local_jacobian_function)
                 {
-                    std::cerr << "The function 'local_jacobian_function' of operator '" << scheme.name() << "' has not been implemented."
-                              << std::endl;
+                    samurai::io::eprint("The function 'local_jacobian_function' of operator '{}' has not been implemented.\n", scheme.name());
                     assert(false && "Undefined 'local_jacobian_function'");
                     exit(EXIT_FAILURE);
                 }
@@ -114,8 +114,8 @@ namespace samurai
             {
                 if (!m_unknown)
                 {
-                    std::cerr << "Undefined unknown for this non-linear system. Please set the unknowns using the instruction '[solver].set_unknown(u);'."
-                              << std::endl;
+                    samurai::io::eprint(
+                        "Undefined unknown for this non-linear system. Please set the unknowns using the instruction '[solver].set_unknown(u);'.\n");
                     assert(false && "Undefined unknown");
                     exit(EXIT_FAILURE);
                 }
@@ -264,7 +264,7 @@ namespace samurai
                 }
 
                 // MatView(B, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                // std::cout << std::endl;
+                // samurai::io::print("\n");
 
                 VecRestoreArrayRead(x, &x_data);
 
@@ -285,14 +285,14 @@ namespace samurai
                     using namespace std::string_literals;
                     const char* reason_text;
                     SNESGetConvergedReasonString(snes, &reason_text);
-                    std::cerr << "Divergence of the non-linear solver ("s + reason_text + ")" << std::endl;
+                    samurai::io::eprint("Divergence of the non-linear solver ({})\n", reason_text);
                     // VecView(b, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                    // std::cout << std::endl;
+                    // samurai::io::print("\n");
                     // assert(check_nan_or_inf(b));
                     assert(false && "Divergence of the solver");
                     exit(EXIT_FAILURE);
                 }
-                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout << std::endl;
+                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); samurai::io::print("\n");
             }
 
           public:
diff --git a/include/samurai/petsc/nonlinear_solver.hpp b/include/samurai/petsc/nonlinear_solver.hpp
index f8dd96929..5c2f31479 100644
--- a/include/samurai/petsc/nonlinear_solver.hpp
+++ b/include/samurai/petsc/nonlinear_solver.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../print.hpp"
 #include "fv/cell_based_scheme_assembly.hpp"
 #include "fv/flux_based_scheme_assembly.hpp"
 #include "fv/operator_sum_assembly.hpp"
@@ -133,8 +134,8 @@ namespace samurai
 
                 if (assembly().undefined_unknown())
                 {
-                    std::cerr << "Undefined unknown(s) for this non-linear system. Please set the unknowns using the instruction '[solver].set_unknown(u);' or '[solver].set_unknowns(u1, u2...);'."
-                              << std::endl;
+                    samurai::io::eprint(
+                        "Undefined unknown(s) for this non-linear system. Please set the unknowns using the instruction '[solver].set_unknown(u);' or '[solver].set_unknowns(u1, u2...).\n");
                     assert(false && "Undefined unknown(s)");
                     exit(EXIT_FAILURE);
                 }
@@ -220,7 +221,7 @@ namespace samurai
                 }
 
                 // MatView(B, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                // std::cout << std::endl;
+                // samurai::io::print("\n");
 
                 // Put back the real unknown: we need its B.C. for the evaluation of the non-linear function
                 assembly.set_unknown(*real_system_unknown);
@@ -242,7 +243,7 @@ namespace samurai
                 // assembly().set_0_for_useless_ghosts(b);
 
                 // VecView(b, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                // std::cout << std::endl;
+                // samurai::io::print("\n");
                 // assert(check_nan_or_inf(b));
             }
 
@@ -266,14 +267,14 @@ namespace samurai
                     using namespace std::string_literals;
                     const char* reason_text;
                     SNESGetConvergedReasonString(m_snes, &reason_text);
-                    std::cerr << "Divergence of the non-linear solver ("s + reason_text + ")" << std::endl;
+                    samurai::io::eprint("Divergence of the non-linear solver ({})\n", reason_text);
                     // VecView(b, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                    // std::cout << std::endl;
+                    // samurai::io::print("\n");
                     // assert(check_nan_or_inf(b));
                     assert(false && "Divergence of the solver");
                     exit(EXIT_FAILURE);
                 }
-                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); std::cout << std::endl;
+                // VecView(x, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF)); samurai::io::print("\n");
             }
 
           public:
diff --git a/include/samurai/petsc/utils.hpp b/include/samurai/petsc/utils.hpp
index 704c42b76..fc7a37911 100644
--- a/include/samurai/petsc/utils.hpp
+++ b/include/samurai/petsc/utils.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../print.hpp"
 #include <petsc.h>
 #include <xtensor/xfixed.hpp>
 
@@ -160,7 +161,7 @@ namespace samurai
                 {
                     is_nan_or_inf = true;
                     VecView(v, PETSC_VIEWER_STDOUT_(PETSC_COMM_SELF));
-                    std::cout << std::endl;
+                    samurai::io::print(samurai::io::root, "\n");
                     break;
                 }
             }
diff --git a/include/samurai/petsc/zero_block_assembly.hpp b/include/samurai/petsc/zero_block_assembly.hpp
index e0070026e..663ca0b24 100644
--- a/include/samurai/petsc/zero_block_assembly.hpp
+++ b/include/samurai/petsc/zero_block_assembly.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../print.hpp"
 #include "matrix_assembly.hpp"
 
 namespace samurai
@@ -18,7 +19,7 @@ namespace samurai
             {
                 if (value != 0)
                 {
-                    std::cerr << "Unimplemented Assembly(" << value << ")" << std::endl;
+                    samurai::io::eprint("Unimplemented Assembly({})\n", value);
                     exit(EXIT_FAILURE);
                 }
                 this->fit_block_dimensions(true);
diff --git a/include/samurai/print.hpp b/include/samurai/print.hpp
new file mode 100644
index 000000000..82cbf029f
--- /dev/null
+++ b/include/samurai/print.hpp
@@ -0,0 +1,158 @@
+// Copyright 2018-2025 the samurai's authors
+// SPDX-License-Identifier:  BSD-3-Clause
+
+#pragma once
+
+#include <cstdio>
+#include <fmt/format.h>
+#include <fmt/ostream.h>
+#include <type_traits>
+#include <utility>
+
+#ifdef SAMURAI_WITH_MPI
+#include <mpi.h>
+#endif
+
+#include <samurai/arguments.hpp>
+
+namespace samurai
+{
+    namespace io
+    {
+        // Tags & helpers
+        struct all_t
+        {
+        };
+
+        struct root_t
+        {
+        };
+
+        struct rank_t
+        {
+            int value;
+        };
+
+        inline constexpr all_t all{};
+        inline constexpr root_t root{};
+
+        inline constexpr rank_t rank(int v)
+        {
+            return {v};
+        }
+
+        inline int current_rank()
+        {
+#ifdef SAMURAI_WITH_MPI
+            int r           = 0;
+            int initialized = 0;
+            int finalized   = 0;
+            MPI_Initialized(&initialized);
+            if (!initialized)
+            {
+                return 0;
+            }
+            MPI_Finalized(&finalized);
+            if (finalized)
+            {
+                return 0;
+            }
+            MPI_Comm_rank(MPI_COMM_WORLD, &r);
+            return r;
+#else
+            return 0;
+#endif
+        }
+
+        namespace detail
+        {
+            inline bool allow_default_print()
+            {
+#ifdef SAMURAI_WITH_MPI
+                return !(samurai::args::print_root_only && current_rank() != 0);
+#else
+                return true;
+#endif
+            }
+
+            // Default (non-scoped) printing helper
+            template <class... Args>
+            inline void do_print(FILE* s, fmt::format_string<Args...> f, Args&&... a)
+            {
+                if (allow_default_print())
+                {
+                    fmt::print(s, f, std::forward<Args>(a)...);
+                }
+            }
+
+            inline bool should_print(all_t)
+            {
+                return true;
+            }
+
+            inline bool should_print(root_t)
+            {
+#ifdef SAMURAI_WITH_MPI
+                return current_rank() == 0;
+#else
+                return true;
+#endif
+            }
+
+            inline bool should_print(rank_t r)
+            {
+#ifdef SAMURAI_WITH_MPI
+                return current_rank() == r.value;
+#else
+                return r.value == 0;
+#endif
+            }
+
+            // Scoped printing helper
+            template <class Scope, class... Args>
+            inline void do_print(FILE* s, Scope sc, fmt::format_string<Args...> f, Args&&... a)
+            {
+                if (should_print(sc))
+                {
+                    fmt::print(s, f, std::forward<Args>(a)...);
+                }
+            }
+        }
+
+        // stdout
+        template <class... Args>
+        inline void print(fmt::format_string<Args...> f, Args&&... a)
+        {
+            detail::do_print(stdout, f, std::forward<Args>(a)...);
+        }
+
+        // Constrain scoped overloads to known scope tag types to avoid
+        // ambiguity with string literals and other types.
+        template <class Scope,
+                  class... Args,
+                  std::enable_if_t<std::is_same_v<std::decay_t<Scope>, all_t> || std::is_same_v<std::decay_t<Scope>, root_t>
+                                       || std::is_same_v<std::decay_t<Scope>, rank_t>,
+                                   int> = 0>
+        inline void print(Scope sc, fmt::format_string<Args...> f, Args&&... a)
+        {
+            detail::do_print(stdout, sc, f, std::forward<Args>(a)...);
+        }
+
+        // stderr
+        template <class... Args>
+        inline void eprint(fmt::format_string<Args...> f, Args&&... a)
+        {
+            detail::do_print(stderr, f, std::forward<Args>(a)...);
+        }
+
+        template <class Scope,
+                  class... Args,
+                  std::enable_if_t<std::is_same_v<std::decay_t<Scope>, all_t> || std::is_same_v<std::decay_t<Scope>, root_t>
+                                       || std::is_same_v<std::decay_t<Scope>, rank_t>,
+                                   int> = 0>
+        inline void eprint(Scope sc, fmt::format_string<Args...> f, Args&&... a)
+        {
+            detail::do_print(stderr, sc, f, std::forward<Args>(a)...);
+        }
+    } // namespace io
+} // namespace samurai
diff --git a/include/samurai/samurai.hpp b/include/samurai/samurai.hpp
index 181474f7f..1bfbf16b2 100644
--- a/include/samurai/samurai.hpp
+++ b/include/samurai/samurai.hpp
@@ -4,7 +4,6 @@
 
 #ifdef SAMURAI_WITH_MPI
 #include <boost/mpi.hpp>
-#include <fstream>
 namespace mpi = boost::mpi;
 #endif
 #ifdef SAMURAI_WITH_PETSC
@@ -12,11 +11,12 @@ namespace mpi = boost::mpi;
 #endif
 
 #include "arguments.hpp"
+#include "print.hpp"
 #include "timers.hpp"
 
 namespace samurai
 {
-    static CLI::App app;
+    inline CLI::App app;
 
 #ifdef SAMURAI_WITH_PETSC
 #define SAMURAI_PARSE(argc, argv)       \
@@ -67,13 +67,7 @@ namespace samurai
 
 #ifdef SAMURAI_WITH_MPI
         MPI_Init(&argc, &argv);
-        // redirect stdout to /dev/null for all ranks except rank 0
-        mpi::communicator world;
-        if (!args::dont_redirect_output && world.rank() != 0) // cppcheck-suppress knownConditionTrueFalse
-        {
-            static std::ofstream null_stream("/dev/null");
-            std::cout.rdbuf(null_stream.rdbuf());
-        }
+        // No output redirection: fmt wrapper handles output uniformly across ranks
 #endif
         times::timers.start("total runtime");
 
@@ -103,7 +97,7 @@ namespace samurai
         if (args::timers) // cppcheck-suppress knownConditionTrueFalse
         {
             times::timers.stop("total runtime");
-            std::cout << std::endl;
+            samurai::io::print(samurai::io::root, "\n");
             times::timers.print();
         }
 #ifdef SAMURAI_WITH_MPI
diff --git a/include/samurai/schemes/block_operator.hpp b/include/samurai/schemes/block_operator.hpp
index 48e08b557..7588fe1ae 100644
--- a/include/samurai/schemes/block_operator.hpp
+++ b/include/samurai/schemes/block_operator.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../print.hpp"
 #include "../utils.hpp"
 
 namespace samurai
@@ -94,8 +95,11 @@ namespace samurai
                                          using op_field_t = typename operator_t::field_t;
                                          if constexpr (!std::is_same_v<std::decay_t<decltype(u)>, op_field_t>)
                                          {
-                                             std::cerr << "unknown " << i << " is not compatible with the scheme (" << row << ", " << col
-                                                       << ") (named '" << op.name() << "')" << std::endl;
+                                             samurai::io::eprint("unknown {} is not compatible with the scheme ({}, {}) (named '{}')\n",
+                                                                 i,
+                                                                 row,
+                                                                 col,
+                                                                 op.name());
                                              assert(false);
                                              exit(EXIT_FAILURE);
                                          }
diff --git a/include/samurai/schemes/fv/FV_scheme.hpp b/include/samurai/schemes/fv/FV_scheme.hpp
index cdfa0d2e3..109c0cda3 100644
--- a/include/samurai/schemes/fv/FV_scheme.hpp
+++ b/include/samurai/schemes/fv/FV_scheme.hpp
@@ -3,6 +3,7 @@
 #include "../../boundary.hpp"
 #include "../../concepts.hpp"
 #include "../../field.hpp"
+#include "../../print.hpp"
 #include "../../static_algorithm.hpp"
 #include "../../timers.hpp"
 #include "utils.hpp"
@@ -181,7 +182,7 @@ namespace samurai
             {
                 if (input_field.mesh().is_periodic())
                 {
-                    std::cerr << "Error: apply_directional_bc() not implemented for non-box domains with periodic directions." << std::endl;
+                    samurai::io::eprint("Error: apply_directional_bc() not implemented for non-box domains with periodic directions.\n");
                     assert(false);
                     return;
                 }
@@ -320,7 +321,7 @@ namespace samurai
                     return m_dirichlet_config[d];
                 }
             }
-            std::cerr << "No Dirichlet config found for direction " << direction << std::endl;
+            samurai::io::eprint("No Dirichlet config found for direction {}\n", fmt::streamed(direction));
             assert(false);
             return m_dirichlet_config[0];
         }
@@ -373,7 +374,7 @@ namespace samurai
                     return m_neumann_config[d];
                 }
             }
-            std::cerr << "No Neumann config found for direction " << direction << std::endl;
+            samurai::io::eprint("No Neumann config found for direction {}\n", fmt::streamed(direction));
             assert(false);
             return m_neumann_config[0];
         }
diff --git a/include/samurai/schemes/fv/cell_based/cell_based_scheme__nonlin.hpp b/include/samurai/schemes/fv/cell_based/cell_based_scheme__nonlin.hpp
index eb75c38b4..d18cb43ad 100644
--- a/include/samurai/schemes/fv/cell_based/cell_based_scheme__nonlin.hpp
+++ b/include/samurai/schemes/fv/cell_based/cell_based_scheme__nonlin.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../../../print.hpp"
 #include "cell_based_scheme.hpp"
 
 namespace samurai
@@ -178,8 +179,8 @@ namespace samurai
         {
             if (!jacobian_function())
             {
-                std::cerr << "The jacobian function of operator '" << this->name() << "' has not been implemented." << std::endl;
-                std::cerr << "Use option -snes_mf or -snes_fd for an automatic computation of the jacobian matrix." << std::endl;
+                samurai::io::eprint("The jacobian function of operator '{}' has not been implemented.\n", this->name());
+                samurai::io::eprint("Use option -snes_mf or -snes_fd for an automatic computation of the jacobian matrix.\n");
                 exit(EXIT_FAILURE);
             }
 
diff --git a/include/samurai/schemes/fv/explicit_FV_scheme.hpp b/include/samurai/schemes/fv/explicit_FV_scheme.hpp
index 5b6b096f3..f1bc38caf 100644
--- a/include/samurai/schemes/fv/explicit_FV_scheme.hpp
+++ b/include/samurai/schemes/fv/explicit_FV_scheme.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../../print.hpp"
 #include "../explicit_scheme.hpp"
 #include "FV_scheme.hpp"
 
@@ -83,7 +84,7 @@ namespace samurai
 
         virtual void apply(std::size_t /* d */, output_field_t& /* output_field */, input_field_t& /* input_field */)
         {
-            std::cerr << "The scheme '" << scheme().name() << "' cannot be applied by direction." << std::endl;
+            samurai::io::eprint("The scheme '{}' cannot be applied by direction.\n", scheme().name());
             assert(false);
             exit(EXIT_FAILURE);
         }
diff --git a/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_het.hpp b/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_het.hpp
index 1a120fe49..e4ad4ffb0 100644
--- a/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_het.hpp
+++ b/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_het.hpp
@@ -1,5 +1,6 @@
 #pragma once
 // #include "../../../petsc/fv/flux_based_scheme_assembly.hpp"
+#include "../../../print.hpp"
 #include "../explicit_FV_scheme.hpp"
 #include "flux_based_scheme__lin_het.hpp"
 
@@ -63,8 +64,8 @@ namespace samurai
 #ifdef SAMURAI_CHECK_NAN
                                 if (std::isnan(field_value(input_field, comput_cells[c], field_j)))
                                 {
-                                    std::cerr << "NaN detected when computing the flux on the interior interfaces: " << comput_cells[c]
-                                              << std::endl;
+                                    samurai::io::eprint("NaN detected when computing the flux on the interior interfaces: {}\n",
+                                                        fmt::streamed(comput_cells[c]));
                                     assert(false);
                                 }
 #endif
@@ -96,8 +97,8 @@ namespace samurai
 #ifdef SAMURAI_CHECK_NAN
                                     if (std::isnan(field_value(input_field, comput_cells[c], field_j)))
                                     {
-                                        std::cerr << "NaN detected when computing the flux on the boundary interfaces: " << comput_cells[c]
-                                                  << std::endl;
+                                        samurai::io::eprint("NaN detected when computing the flux on the boundary interfaces: {}\n",
+                                                            fmt::streamed(comput_cells[c]));
                                         assert(false);
                                     }
 #endif
diff --git a/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_hom.hpp b/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_hom.hpp
index ecc2cdddc..d8845fc7d 100644
--- a/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_hom.hpp
+++ b/include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_hom.hpp
@@ -1,5 +1,6 @@
 #pragma once
 // #include "../../../petsc/fv/flux_based_scheme_assembly.hpp"
+#include "../../../print.hpp"
 #include "../explicit_FV_scheme.hpp"
 #include "flux_based_scheme__lin_hom.hpp"
 
@@ -287,9 +288,8 @@ namespace samurai
 #ifdef SAMURAI_CHECK_NAN
                                     if (std::isnan(field_value(input_field, stencil.cells()[c], field_j)))
                                     {
-                                        std::cerr
-                                            << "NaN detected when computing the flux on the boundary interfaces: " << stencil.cells()[c]
-                                            << std::endl;
+                                        samurai::io::eprint("NaN detected when computing the flux on the boundary interfaces: {}\n",
+                                                            fmt::streamed(stencil.cells()[c]));
                                         assert(false);
                                     }
 #endif
diff --git a/include/samurai/schemes/fv/flux_based/flux_based_scheme.hpp b/include/samurai/schemes/fv/flux_based/flux_based_scheme.hpp
index 5d41b8fe1..5d903c9c2 100644
--- a/include/samurai/schemes/fv/flux_based/flux_based_scheme.hpp
+++ b/include/samurai/schemes/fv/flux_based/flux_based_scheme.hpp
@@ -1,6 +1,7 @@
 #pragma once
 #include "../../../arguments.hpp"
 #include "../../../interface.hpp"
+#include "../../../print.hpp"
 #include "../../../reconstruction.hpp"
 #include "../../explicit_scheme.hpp"
 #include "../FV_scheme.hpp"
@@ -24,8 +25,9 @@ namespace samurai
         // cppcheck-suppress knownConditionTrueFalse
         if (args::finer_level_flux != 0 && cfg::dim > 1 && cfg::stencil_size > 4 && !args::refine_boundary)
         {
-            std::cout << "Warning: for stencils larger than 4, computing fluxes at max_level may cause issues close to the boundary."
-                      << std::endl;
+            samurai::io::print(
+                samurai::io::root,
+                "Warning: for stencils larger than 4, computing fluxes at max_level may cause issues close to the boundary.\n");
         }
 
         return FluxBasedScheme<cfg, bdry_cfg>(flux_definition);
diff --git a/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp b/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp
index c8c5eced7..b5a78f236 100644
--- a/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp
+++ b/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "../../../print.hpp"
 #include "flux_based_scheme.hpp"
 
 namespace samurai
@@ -103,8 +104,9 @@ namespace samurai
         {
             if (enable && dim > 1 && stencil_size > 4 && !args::refine_boundary) // cppcheck-suppress knownConditionTrueFalse
             {
-                std::cout << "Warning: for stencils larger than 4, computing fluxes at max_level may cause issues close to the boundary."
-                          << std::endl;
+                samurai::io::print(
+                    samurai::io::root,
+                    "Warning: for stencils larger than 4, computing fluxes at max_level may cause issues close to the boundary.\n");
             }
             m_finer_level_flux = enable ? -1 : 0;
         }
@@ -595,8 +597,8 @@ namespace samurai
 
                 if (!jacobian_function)
                 {
-                    std::cerr << "The jacobian function of operator '" << this->name() << "' has not been implemented." << std::endl;
-                    std::cerr << "Use option -snes_mf or -snes_fd for an automatic computation of the jacobian matrix." << std::endl;
+                    samurai::io::eprint("The jacobian function of operator '{}' has not been implemented.\n", this->name());
+                    samurai::io::eprint("Use option -snes_mf or -snes_fd for an automatic computation of the jacobian matrix.\n");
                     exit(EXIT_FAILURE);
                 }
 
diff --git a/include/samurai/stencil.hpp b/include/samurai/stencil.hpp
index c2a2aca46..466564cc8 100644
--- a/include/samurai/stencil.hpp
+++ b/include/samurai/stencil.hpp
@@ -1,5 +1,6 @@
 #pragma once
 #include "indices.hpp"
+#include "print.hpp"
 #include "static_algorithm.hpp"
 
 namespace samurai
@@ -591,7 +592,7 @@ namespace samurai
 #ifndef NDEBUG
             if (origin_cell.index > 0 && static_cast<std::size_t>(origin_cell.index) > m_mesh.nb_cells()) // nb_cells() is very costly
             {
-                std::cout << "Cell not found in the mesh: " << origin_cell << std::endl;
+                samurai::io::eprint("Cell not found in the mesh: {}\n", fmt::streamed(origin_cell));
                 assert(false);
             }
 #endif
@@ -628,7 +629,9 @@ namespace samurai
 #ifndef NDEBUG
                         if (cell.index > 0 && static_cast<std::size_t>(cell.index) > m_mesh.nb_cells()) // nb_cells() is very costly
                         {
-                            std::cout << "Non-existing neighbour for " << origin_cell << " in the direction " << dir << std::endl;
+                            samurai::io::eprint("Non-existing neighbour for {} in the direction {}\n",
+                                                fmt::streamed(origin_cell),
+                                                fmt::streamed(dir));
                             // save(fs::current_path(), "mesh_error", {true, true}, m_mesh);
                             exit(1);
                         }
diff --git a/tests/test_operator_set.cpp b/tests/test_operator_set.cpp
index e14b026cf..078488e78 100644
--- a/tests/test_operator_set.cpp
+++ b/tests/test_operator_set.cpp
@@ -18,7 +18,7 @@ namespace samurai
 
     //     LevelCellArray<1> lca{lcl};
 
-    //     std::cout << lca << "\n\n";
+    //     samurai::io::print("{}\n\n", fmt::streamed(lca));
     // }
 
     template <std::size_t size>